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ABSTRACT 


Improved control of gas turbine propulsion plants could 
offer the Navy increased economic, maintenance, and tactical 
benefits. This thesis provides methods of steady state and 
dynamic computer modelling, as well as two non-proprietary 
control design methods. The classical proportional integral 
(PI) regulator design method and the modern linear quadratic 
regulator (LQR) controller design method are used to demon- 
strate a base for Navy redesign of existing gas turbine 
propulsion plant control systems. A comparison of the PI 
and LOR designs is conducted. In addition, a real or near- 
real time dynamic computer simulation is presented that has 
immediate application in the areas of model-based control 


and plant health monitoring. 
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I. INSRODUCTION 


Many modern surface Naval marine propulsion plants are a 
combination of gas turbine engines with controllable 
reversible pitch propellers. This presents the problem of 
matching the engine RPM to the most efficient pitch which 
has been accomplished through the use of an integrated 
throttle control (ITC). Figure 1 shows a block diagram of a 
typical ITC control scheme. 

The implementation technology for Figure 1 is well over 
20 years old and its limitations are now well defined. 
Today, technology exists that will allow the antiquated 
analog mechanisms and current computerized systems to be 
replaced by smaller more reliable digital controls and 
hardware. This approach suggests that the following 
benefits could be realized: 

(1) Reduction of maintenance "nightmares" that develop 
due to the intricacy and number of small parts in 


components such as mechanical fuel governors; 


(2) More reliable and compact circuitry would modify 
present hardware such as the Free Standing Electronic 


Enclosure, propulsion and electrical control 
consoles, and current engine health monitoring 
equipment; 


(3) Advances in the ability to model and simulate gas 
turbine performance would allow plant performance to 
be significantly improved, thereby increasing plant 
efficiency and translating to lower operating costs; 


(4) New techniques in engine health monitoring and 
analysis provide essential real time data on plant 
performance to the operators, allowing better and 
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more rapid evaluation and response to a potential or 
actual engineering casualty; 


(5) More compatibility between control systems could be 
achieved, thereby reducing the number of different 
repair parts that must be stocked in the naval supply 
system. More commonality would also streamline the 
training process of personnel responsible for 
maintaining and operating the systems; 

(6) Inherent flexibility through reprogramming of 
computer-based controls would pave the way for future 
developments. 

A. CONTROLLER BACKGROUND 

During the late seventies and early eighties the marine 
gas turbine industry hotly debated the pros and cons of 
analog vs. digital control to implement integrated throttle 
control [{Ref. 1]. The advocates of analog control were of 
the opinion that this technology was reliable and could 
perform all necessary calculations required for effective 
mpant control. It was felt that little would be gained in 
the way of reduction of component count or system reliabili- 
ty through digital systems. This thought process led Rolls 
Royce to choose analog systems for warship controls, and led 
General Electric to a similar conclusion for the main fuel 
control on the LM-2500. 

The digital advocates on the other hand, had the 
foresight to realize that advances in technology would be 
more easily implemented in a digital base, and that 
reliability would indeed be as good, if not better than, 


analog systems. With the advent of the microprocessor, the 


component count can indeed be reduced with a carefully 


executed design process. This was demonstrated by the 
aviation community first on the F100 engine (Ref. 2}. A 
natural progression would be for the marine gas turbine 
conmunity to feilow sean. It must be realized that some 
analog fuel system control components will probably always 
be required, particularly in the sensing and actuation 
areas. 

Perhaps the most compelling reason today to convert to 
digital control is the advent of intelligent control. Tem 
this approach, it is possible to control a large quantity of 
measured and unmeasured variables with a limited amount of 
operator intervention to meet the dynamic needs of the 
operator. 

Typically, a good control design approach consists of 
eleven steps. These steps contain three "feedback loops" 
which provide the means for modification or improvement 
should the deSigner desire. This control deSign approach is 
as follows: 

(1) Specifications for control design; 

(2) Evaluation of plant function; 

(3) Plant mathematical modeling; 

(4) Plant model validation--open loop simulation; 
(5) Selection of control strategy; 

(6) Selection of actuators, sensors; 

(7) Dynamic modeling of actuators, sensors; 


(8) Selection of controller action; 


(9) Theoretical controller design; 
(10) Controller validation--closed loop simulation; 
(11) Prototype. 

The design feedback loops exist between steps 4 and 3, 
between steps 10 and 8, and between steps 10 and “5. 
Utilization of computer-aided design techniques in the 
design, validation, and optimization of control schemes 
provides an efficient and economical method for selecting 
the most suitable candidate for hardware development. 
Prudent selection of designs is essential considering the 
complexity and large capital expenditure incurred as the 
design progresses from the conceptual phase to its final 
form. Inherent in this approach is the need for evaluation 
and modelling of gas turbine performance (step 3). Conse- 
quently, while this thesis is dealing with marine gas 
turbines, much early work was done in the area of aviation 
gas turbine modeling and control. We begin with a review of 
these efforts. Chapter II is a review of previous recent 
work in gas turbine modelling and control. Chapter III is a 
review of work previously performed at the Naval Postgradu- 
ate School. Chapters IV and V detail the steady state and 
dynamic simulations of this work. Chapter VI contains the 
design methods for a classical proportional integral (PI) 
regulator and a modern linear quadratic (LOR) regulator. 


Also in Chapter VI the controller designs are compared for 


operation in a simulated sea _ state. Conclusions and 


recommendations make up Chapter VII. 


II. PREVIOUS GAS TURBINE MODELLING AND CONTROL 


A. EARLY COMPUTER MODELS 
Gas turbines in use today for marine propulsion are for 


the most part derivatives of aviation gas turbine engines 


that have been "marinized" for use at sea. As one would 
expect, several computer simulations were developed to 
evaluate and predict system performance. The early 


simulations were developed by the aviation industry anda 
provided a substantial data base for development of more 
advanced computer models. A short summary of some of the 
major early aircraft simulations is given below [Ref. 3]: 


(1) SMOTE--Developed in 1967 by the Turbine Engine 
Division of the U.S. Air Force Propulsion Laboratory 
(AFAPL), Wright-Patterson AFB, Ohio. It us capable 
of calculating steady-state design and off design 
performance of a two-spool turbofan engine. 


(2) GENENG--Developed in 1972 by NASA's Lewis Research 
Center (LeRC), Cleveland, Ohio. TtSs purpose is to 
improve the versatility of SMOTE. Steady-state 
design and off-design performance of one- and two- 
spool turbojets can be calculated as well as the two- 
spool turbofan. 


(3) GENENG JII--Derivative of GENENG, it calculates 
steady-state performance of two- or three-spool 
turbofan engines with as many as three nozzles. 


(4) NEPCOMP--Developed in 1974 by the Naval Air 
Development Center (NADC), Warminister, Pennsylvania. 
The flexibility inherent to NEPCOMP allows’ for 
calculation of steady-state performance of gas- 
turbine engines with multiple spools, including 
turbojets, turbofans, turboshafts, and ramjets. 


(5) DYNGEN--Developed in 1975 by LeRC, it combined the 
capabilities of GENENG and GENENG II for calculating 


steady-state performance of gas turbine engines with 
multiple spools. The additional capability of 
calculating transient performance was also added. 

(6) NNEP--Jointly developed in 1975 by NASA, LeRC, and 
NADC. This computer code is able to simulate steady- 
state design and off design performance of almost any 
conceivable gas turbine engine simulation. 

As can be seen above, the majority of the early work was 
devoted to steady-state simulations. A major shortfall was 
a lack of dynamic simulation capability. At this point it 
is prudent to shift the emphasis from the work performed by 
the aviation industry and concentrate on the contributions 


made in the marine gas turbine industry in the area of 


dynamic simulation. 


Be DYNAMIC COMPUTER MODELS FOR MARINE ENGINE SIMULATION 

Engineers at David W. Taylor developed equations to 
mathematically model various engine components for a 
"building block" approach to modelling [Ref. 3]. Once these 
were established, a system of common component interface 
locations was defined and the locations were numbered. 


Equations were then developed for the numbered major gas 


turbine components, including compressors, burners, 
turbines, and engine load. Dynamic equations were then 
developed to describe speed, power balances, mass 
accumulation, and energy accumulation.1t An iterative 


approach was then utilized to balance the performance 


linformation used for this portion of the discussion 
only relates to a simulation of a single spool engine 
ConlErcirae tom. 


characteristics of the various engine components. A Newton- 
Raphson technique was used to achieve convergence. The 
results of the simulations yielded good comparisons between 
the manufacturer's simulation and the existing experimental 
data. 

Beginning in the early seventies, the U.S. Navy 
initiated The Gas Turbine Ship Propulsion Control Systems 
Research and Development Program [Ref. 4]. The Navy chose 
Propulsion Dynamics, Incorporated to conduct the program 
which was designed to develop a machinery dynamics’ and 
control system data base. The program involved computer 
Simulations of total propulsion systems, which were 
validated by shipboard and model testing. The program 
continued into the eighties and was still generating 
technical papers as recently as 1986. The program was 
successful in developing a theoretical design base for gas 
turbine propulsion systems. Major conclusions were drawn in 
the following areas [{Ref. 5]: 

(1) Propulsion systems cycling; 

(2) Propeller speed governing; 

(3) Gas generator power governing; 

(4) Combined Power and Speed Governing. 

Based on data obtained during the program, ae ship 
propulsion control system was devised for use in computer 
Simulations. The control system was of the classical 


integral variety, whose gains were fixed via a "cut and try" 


method. Linear controller gains were obtained for various 
wave conditions and engine speeds, then tabulated and 
compared. In a current application the gain is fine tuned 
via the "sea state adjust" control found on the propulsion 
control consoles aboard DD-963 class destroyers to account 
for variations in the load and non-linear propulsion plant. 
Figure 1 shows a block diagram of the ship propulsion 
control system used. Simulations of this approach tended to 
give good results when compared to model and ship generated 
data\"[Rebk. Si): However, the approach generated some 
interesting observations regarding a gas turbine engineering 
plant's response to seaway- and maneuver-induced unsteady 
loading, which are indeed confirmed by the experience of the 
author who served as Main Propulsion Assistant aboard a DD- 
963 class destroyer. In high seas, gas turbine plants 
experience a good deal of engine/propeller cycling due to 
constant changes in propeller loading as the ship moves 
through the water. A ship configured with two propulsion 
shafts experiences a good deal of propeller load variation 
during turns, particularly during high speed turns. 
Naturally, these conditions cause numerous changes in engine 
speed, resulting in engine wear and potential overspeeding 
of the engine gas generator should the propulsion load be 
lost for some reason. It should be noted at this point that 
these two phenomena can be thought of as "disturbances" to 


the plant. 
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Returning to general control development, modern control 
theory provided the next logical step in controller design. 
In this work, state space techniques applied to gas turbines 
have yielded positive results. Such state variable methods 
allow the control system designer to gain an understanding 
of the inherent input cross-channel coupling dynamic 
characteristics of the system and to take advantage of 
coupling which exists between input and output variables. 

In the late seventies students and faculty at the Naval 
Postgraduate school applied state space techniques to a 
linearized model of an FFG-7 ship propulsion system [Ref. 
6]. Dynamic propulsion system equations were developed for 
the FFG-7 and then linearized, the appropriate matrices 
developed, and the dynamic simulations conducted. The 
results demonstrated that the linear model described the 
system behavior reasonably well. 

Another mathematical model was developed at Tsinghua 
University, Beijing, China in the mid-eighties [Ref. 7]. A 
three shaft marine gas turbine was modelled and simulated 


using state space techniques, and two different numerical 


methods were used to obtain convergence. The convergence 
methods used were: (1) the varying coefficient method, and 
(2) the small deviation method. The difference in methods 


lies in the fact that only small system perturbations can be 
considered in the latter, while large perturbations can be 


considered in the former. In the first method the initial 


ea 


point of linearization lies in the unsteady regime. The 
real beauty of the varying coefficient method is that 
transients under large perturbations can be obtained with 
sufficient accuracy using linearized equations. Results 
from the two simulation techniques were compared and the 


varying coefficient method was deemed more accurate. 


C. RECENT CONTROL DESIGN TECHNIQUES 

There are numerous methods by which one can design a 
modern controller for an automatic system. When a state 
space approach is taken to design, there are basically two 
ways to approach the task: (1) The Pole Placement method, 
and (2) The Linear Quadratic Regulator technique (LQR). 
The Pole Placement method requires that the location of the 
desired system closed loop poles be known. Since the 
optimum closed loop poles of a system may not be known 
during design, the LQR method is often a better choice. The 
LOR method optimizes the design of the controller, based on 
the inputs of various matrices and a cost function. The LOR 
controller often requires an observer to calculate the 
states, it then calculates the error between actual and 
desired states and computes the gains such that stability is 
guaranteed and the integrated error minimized. (The theory 
of this approach will be reviewed in more detail in the 
following chapters). 

Kidd, Munro, and Winterbone examined the potential of a 


digital control scheme designed using LOR state space 
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techniques [{Ref. 8]. The plant model was one of a two- 
shaft, two-turbine vessel with a combination of a sprint and 
a cruise turbine on each shaft coupled to a controllable 
reversible pitch propeller via a reduction gear. The 
simulations were performed using a FORTRAN IV digital, non- 
linear, dynamic computer simulation which included steady 
state data for the non-linear propeller and thrust 
characteristics. A digital controller was developed using 
state space techniques, eventually culminating in a gain- 
scheduled multivariable controller which was’ constructed 
from a selection of linear compensator designs. For 
comparison purposes a conventional control system was 
designed as a yardstick by which to meaSure the digital 
control system. Both controllers were then implemented in 
the non-linear ship simulation model. The responses of the 
two controllers were compared for several maneuvers and the 
multivariable controller demonstrated a much faster speed of 
response and less overshoot on propeller-shaft torque 
SMe put . The multivariable controller constrained the 
propeller well within safe and acceptable operating limits. 
The improvements in response of the propulsion plant 
improved the ship speed response which resulted in ship 
acceleration and stopping time improvements, i.e ship 
maneuverability improvements. 

LQR controllers have also been designed for the F-401 


and F100 aerospace turbofan engines. Figure 2 is a block 


3 


diagram of the F100 control model [(Ref. 2]. Similar 
research was done to apply LOR techniques to the design of a 
power turbine governor for a turboshaft engine driving a 
helicopter rotor blade [Ref. 9]. In that work, a GE-700 
turboshaft engine was modelled using state space methods and 
was mathematically coupled to a linear lumped capacitance 
model of an articulated rotor blade. The two were then 
combined into an overall system matrix and simulated; the 
results were compared to a conventional governor's 
performance. The performance was increased in the areas of 
time response and overshoot in power turbine speed. These 
results seem to parallel the results obtained by Kidd, 


Munro, and Winterbone, but for a different application. 
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III. PREVIOUS WORK AT NPS 


The plant being considered here is a Boeing Model 502- 


6A 175 horsepower gas turbine connected to a Clayton 17-300 


water brake dynamometer as shown in Figure 3. The gas 
turbine is divided into two separate sections, a gas 
generator section and a power output section. The gas 


generator is composed of a single entry, single stage 
centrifugal compressor connected to a single stage axial 
flow high pressure turbine (HPT). Two cross-connected 
through-flow type combustion chambers provide an aerodynamic 
coupling between the turbine and compressor. An accessory 
drive section is geared off the gas generator shaft, and 
contains the electric starter, tachometer generator, fuel 
pump, governor, and lube oil pump. The power output section 
consists of an axial flow free power turbine (FPT), reduc- 
tion gearing, and output shaft. The gas generator and power 
output sections are connected aerodynamically. 


Previous work by Johnson [{Ref. 10] (hardware design and 


implementation), Herda [Ref. 11] (computer modelling and 
simulation), and Miller {[{Ref. 12] (model testing and 
modification), has provided the starting point for the 


present work. 
The state space model previously developed has been 


slightly modified, but remains essentially the same with the 
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exception of one of the plant inputs. The model for the 


present work applies the state space linearization: 
x = Akx + Btu (1) 


where 


x = the state vector, 
u = the input vector, 
A = the state coefficient matrix, 
B = The input coefficient matrix. 


The states are defined as gas generator speed (NG), 
power turbine/dynamometer speed (NS), and mechanical energy 
resulting from fuel combustion (E). The plant inputs are 
fuel flow rate (MF) and dynamometer torque (QD). 

Perturbations which are the basis for the linearizations 


are accomplished by using the equation: 


Kame ie a (2) 


where 


Xo — thee sy aes 


m* 
Il 


§X = the perturbation from the initial value, 


X = the current value. 
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All plant state space variables are represented in this 
manner. Employing the perturbational variables, the state 


space equation becomes: 


ng ng mf 
ns = A ns ate B qd 
e e 


a 


The elements of the "A" and "B" matrices are calculated 
using Taylor serieS expansions on each plant component, 
retaining only first order terms. So, the elements of the 


state space "A" and "B" matrices can be written symbolically 


as: 
all = 93Ng/dng al2 = 3Ng/dns al3 = dNg/ce 
a21 = 9Ns/ dang a22 = dNs/ons a23 = dNS/de 
a31 = 9E/9ng a32 = dE/ons a33. = 0 5E/ de 
b11 = 3Ng/amf b12 = 3Ng/dqd 
b21 = 3Ns/omf b22 = 3Ns/dqd 
b31 = 3E/omf b32 = 3E/3qd 
Herda developed both steady state and dynamic computer 
Simulations to describe the behavior of the plant. The 


dynamic equations were derived from quadratic curve fits of 
steady state data collected during operation of the gas 
turbine/dynamometer. Uncorrected variables were used, 
primarily because the conditions in the test cell remained 
near standard conditions at all times. The modifications to 


the computer code to correct for temperature and pressure 
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could easily be added after successful concept validation. 
The cause and effect multiport model used for that work is 
depicted by Figure 4. It was initially proposed by Johnson, 
then expanded by Herda. 

It is apparent from the multiport model that the 
variable coupling the gas generator and power output 
sections was the pressure P4. P4 is both the high pressure 
turbine exhaust pressure and the free power turbine inlet 


pressure. Herda's steady state model was developed on the 





premise that for any operating point there was one fuel flow 
rate (MF), one high pressure turbine inlet pressure (P2), 
and one high pressure turbine exhaust/free power turbine 
inlet pressure (P4). Inputs to the model were gas generator 
speed (NG) and dynamometer speed (NS). From these inputs an 
initial fuel guess was computed. Convergence to the correct 
P2 and P4 was then attempted using a modified Newton-Raphson 
algorithm. If the pressures could not be converged, the 
fuel estimate was modified using a golden section technique. 
These iterations continued until either satisfactory 
convergence to specified criteria was obtained (in which 
case a torque comparison between the compressor and high 
pressure turbine was performed) or convergence failed and no 
solution was obtained. If the torque comparison failed to 
meet its convergence criteria, the iterations continued as 
just described. If satisfactory convergences between the 


pressures and torques were obtained, the routine calculated 
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the remaining plant variables and the state space "A" and 
"B" matrices. 

Herda performed limited testing of his steady state 
model and was able to obtain some good results, but he did 
experience numerical instability and failed convergence in 
some instances. Herda's dynamic model was developed using 
mainframe dynamic simulation software. 

Miller's work focused on solving the nhuUmendgean 
instability problem encountered by Herda. A performance 
envelope was developed and additional data obtained for 
analysis. Miller made minor modifications to the model and 
investigated alternative convergence methods with some 
success, but he also encountered numerical instabilities. 

A summary of previous work is as follows: 


(1) The cause-effect multi-port model worked well to 
portray system response; 


(2) Computer algorithms derived from the multi-port model 
provide a method for linearization and state space 
matrix computation based upon steady state experimen- 
tal data; 


(3) Numerical convergence for the steady state model was 
uncertain for portions of the plant operating 
envelope; and, 


(4) Great accuracy was required in the computation of 


steady state plant variables, particularly the 
pressure P4. 


Ze 


IV. STEADY STATE SOLUTION 


The ability to obtain a steady state solution at any 
point in the operating envelope was prerequisite to the 
later dynamic simulation and controller design phases of 
this work. The criteria for an acceptable steady state 
solution was to achieve a torque balance between the 
compressor and high pressure turbine. A zero torque 
aifferential is indicative of gas generator balance and 
steady state operation, while a non-zero torque differential 
1s indicative of gas generator acceleration or deceleration. 

As in Herda's solution method, there was assumed to be a 
distinct MF, P2, and P4 for every steady state torque value, 
and these quantities were required for calculation of the 
remaining plant variables. 

Techniques previously used for converging torque and 
pressure values were slope methods. Several other 
mathematically intense slope convergence algorithms were 
initially investigated in the present work, these also 
proved to be inadequate. To better visualize the behavior 
of the torque and pressure curves being dealt with, a torque 
balance equation was derived using the compressor and high 
pressure turbine equations developed by Herda. The 
resulting equation was a quadratic expression in terms of 


five variables: MF, P2, P4, NG, and NS. Three of these 


23 


(MF, NG, and NS) were known or could be calculated for a 
particular operating point, leaving P2 and P4 as unknowns. 
A grid procedure was employed to solve for P2 and P4 which 
forced two criteria to be met: 

(1) The gas turbine must be torque balanced; 


(2) All component input-output relationships must be 
satisfied. 


The modified multiport diagram of Figure 5 illustrates the 
process. 

Fuel (MF) and a range of potential P2G values (P2 
"guessed") were used as inputs. A corresponding value of 
P4G was calculated from the torque balance quadratic 
equation. An imaginary root check discarded any imaginary 
roots, leaving only the real roots for consideration as 
possible solutions. Roots acquired using the negative 
radical portion of the quadratic equation were defined as 
"low energy" solutions, while roots acquired using the 
positive radical were defined as "high energy" solutions. 
It was decided that should the situation arise where both 
high and low energy solutions existed in P4G, the low energy 
solution would be chosen because it corresponded physically 
to less fuel used for the operating point. Each” pair gem 
guessed pressures (P4G) was then input into calculations to 
check for torque balance. If torque balance was achieved, 
the corresponding values were recorded and the routine 
continued. If the torque balance checks failed, the next 


value of P2G was input and the process repeated. Torque 
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balance caused the computation of P4C, the computed P4 
pressure, which was calculated from a subroutine involving 
component input output relationships with P2G and MF as the 
inputs. Crossing logic was then used to detect points where 
the P4C and P4G curves met (or crossed). A crossing of 
these curves was considered to be a potential operating 
point of the plant. Figure 6 showS an example of P4C and 
P4G curves crossing. 

The results of this procedure fell into one of the 
following categories: 


(1) A solution existed, but outside of the valid P2 
range; 


(2) Multiple crossings existed in the valid P2 range; 
(3) Imaginary solutions existed; 

(4) A combination of 1, 2, and 3; 

(5) Only one solution existed in the valid P2 range; 
(6) No solution existed (no crossings). 

Convergence to a root in category 1 or 3 above could 
result in an incorrect steady state solution. This 
explains some of the numerical instability and inability to 
converge some points in the operating envelope encountered 
in the past. 

The existence of multiple roots presents the rather 
formidable task of consistently extracting the root that 
leads to the correct steady state solution. In the case of 


two roots in the valid P2 range, it was discovered that 
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often the high energy solution provided a better steady 
state answer than the low energy solution. 

Eventually this method was abandoned because of 
inaccuracies in equation coefficients and the large changes 
in the "A" and "B" matrices which occurred for very small 
changes in P4. For this reason a grid search technique was 
adopted. Although computationally intense, the method 
considered all possible combinations of MF, P2, and P4 
within specified ranges, thereby eliminating the problem of 
converging on the first root which occurs in gradient 
methods. 

Two computer algorithms were developed, one to provide 
the MF, P2, and P4 ranges to be searched, and the second to 
converge these three values. Figure 7 is a flowchart for 
the first algorithm which is called the Variable Range 
Determining (VRD) algorithm, a copy of which is included as 
Appendix A. The inputs to the program were gas generator 
speed (NG), power turbine speed (NS), and the number of 
iterations for each of the variables: MF, P2, and P4. This 
number of iterations determined both the size of the search 
increments and the width of search area. The fuel initial 
guess was computed by subroutine NGNSMF to determine the 
starting point for fuel iteration. A copy of NGNSMF and all 
other subroutines used is included as Appendix D. The 


variable initialization section was then entered to set the 
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Figure 7. Flowchart of VRD Algorithm 
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values of the various increments as well as to initialize 
the high and low values. The VRD increments were 
established using the maximum and minimum possible values in 
the normal plant operating envelope for each of the three 
variables being considered. The first fuel guess was 
arbitrarily set 20 pounds less than the value returned by 
NGNSMF to ensure proper coverage of the beginning of the 
fuel range. The fuel iteration loop then incremented the 
fuel and set P2 to the low value of its operating range. 
The P2 loop was then entered, P2G incremented, and 
subroutines called to compute compressor torque (OGim 
compressor discharge temperature (T2), and air mass flow 
rate (MA). The P4 low value was then set, P4 loop entered, 
and P4G incremented. The high pressure turbine torque 
(QHPT) was calculated via a subroutine, and then QC and QHPT 
were compared. If the difference had not changed sign, 
torque convergence had not occurred and P4G was incremented. 
If the difference had changed sign or was equal to zero, 
torque convergence waS assumed and subroutines to calculate 


high pressure turbine discharge temperature (T4) and P4C 


were called. P4G and P4C were then compared for convergence 
in a similar fashion. If convergence was not obtained, P4G 
was incremented and the process continued. If P4 conver- 


gence was obtained the P2 subroutine was called upon to 
compute P2C, then the third and last convergence check was 


made on  P2. Failure to converge incremented P4G. 
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Successful convergence in all three tests resulted in logic 
which identified the variable range for possible 
solution(s). Upon completion of the P4 range, the P2 loop 
was incremented and the process repeated. The MF loop was 
incremented when the range of P2 values was exhausted, and 
the routine continued until the MF range had been traversed. 
The end result was that for every incremented value of MF, 
all combinations of P2G and P4G were examined and compared 
to computed values. The results were then grouped to 
provide the solution variable ranges for the second 
algorithm, which refined the solution. 

The Solution Convergence (SC) algorithm was essentially 
the same as the previous VRD algorithm. A copy of the SC 
program iS included as Appendix B. The high and low values 
for MF, P2, and P4 were specified to be those obtained from 
the VRD algorithm. The increments in the SC algorithm were 
defined as functions of the VRD high and low variables and 
the originally defined VRD increments. The SC increments 
were smaller than the VRD increments, providing a finer grid 
to be searched. The search range for each variable was 
established by starting one VRD increment outside the 
initial value and then incrementing by SC increments. This 
method ensured proper coverage in the vicinity of the 
initial value of the range. Coverage was extended slightly 
past the final value by adding an arbitrary number of 


iterations to the number of "guesses" specified for each 


Sek 


variable during the input phase. Convergence logic was the 
Same as previously discussed. The accuracy of the converged 
solution was determined by the magnitude of the terms DELQ 
(QC - QHPT), DELP2 (P2C - P2G), and DELP4 (P4C - P4G). The 
smaller the "DEL" terms, the more accurate the solution. 
The output of the routine was a list of all converged 
solutions that lay within the specified ranges of MF, P2, 
and P4. 

With the critical convergence criteria met, the final 
portion of the steady state solution process could be 
completed. A third computer routine was developed from work 
performed by Herda and Miller. The Steady State Variable 
and Partial Derivative (SSVPD) Algorithm used the converged 
MF, P2, and P4 values to compute steady state values for air 
mass flow (MA), air-fuel mass flow (MAF), high pressure 
turbine discharge temperature (T2), power turbine discharge 
temperature (T4), free power turbine torque (QFPT), 
dynamometer torque (QD), high pressure turbine torque 
(QHPT), compressor torque (QC), and dynamometer water weight 
(WW). SSVPD calculated the partial derivatives required for 
the necessary linearizations to form the state space "A" and 
"B" matrices. A copy of the SSVPD program is included as 
Appendix C. 

Using this three step process, it was possible to 
converge steady state solutions for all gas generator speeds 


between 22000 RPM and 32000 RPM, and for dynamometer speeds 


Je 


between 500 RPM and 2500 RPM. Gas generator speeds below 
22000 RPM were sporadically convergeable, they were 
unconvergeable below 20000 RPM and above 35000 RPM. Power 
turbine speeds above 2500 RPM were not’ consistently 
convergeable. 

Miller hypothesized that quadratic data fits to data for 
the various components of the model may not be valid 
throughout the operating envelope, and this work seems to 
validate that theory. That is, quadratic fits appear to be 
reasonable in the middle of the operating envelope, but not 
in the low or high portions. 

A Steady State Convergence Map of "A" matrices was 
constructed for the convergeable region of the operating 
envelope, and is shown in Figure 8. Since the states NG 
and NS characterize the plant in state space, these 
variables were chosen as the coordinate axes of the grid. 
For each node of the grid, the list of converged solutions 
from the SC routine was examined and the DELQ, DELP2, and 
DELP4 values compared. Convergence accurate to 0.1 pound of 
fuel and 0.1 psi for both pressures were set aS minimums or 
the solution was eliminated from the list. All remaining 
candidates for each node were subsequently entered into the 
SSVPD algorithm and the results collected. Strip Chart 
recordings of actual plant runs were examined and those runs 
lying in the convergeable region were utilized. MThe start 


and end points of these runs were converged and used to 
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anchor the grid. Once anchored, the remaining grid points 
were selected by comparing the various matrix coefficients 
for trends both horizontally along lines of constant NG and 
vertically along lines of constant NS. The sensitivity to 
convergence accuracies was demonstrated by the wide variance 
in magnitude of the matrix coefficients for a given grid 
point, particularly the Al3 and A23 entries. The Al13 entry 
was extremely sensitive to P4. By establishing the 
horizontal and vertical trends on the grid, it was a 
relatively simple matter to select/adjust a particular grid 
point from the various candidates until the best overall 


result was obtained. 
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V. DYNAMIC SIMULATION 


Due to very low gas turbine inertias, the smoothness of 
the changes in "A" matrix entries was critical to effective 
plant “simuloaewone Consequently, the Steady State 
Convergence Map was first analyzed for horizontal and 
vertical trends, both in magnitude and sign. Overall trends 
were readily apparent with the exception of seven Al2 
entries, five A21 entries, and five A23 entries. Of these 
discrepancies one was an ill fitting data point (the A23 
entry of the matrix at NG = 23150 rpm and NS = 493 rpm was 
of the wrong magnitude) and consequently the entire matrix 
was disregarded, while the remainder were of the correct 
magnitude but of the wrong sign. Table 1 is a summary of 
the matrix coefficients that were of the wrong Sign. 
Possible reasons for these errors include poor data fit at 
low engine and dynamometer speeds, and the unreliability of 
data values recorded at low engine horsepowers as documented 
in the dynamometer technical manual (13). Also, these 
values were observed to be extremely sensitive to P4 which 
had a large convergence tolerance. 

Further examination of the Steady State Convergence Map 
revealed that relatively smooth curves could be generated 
both horizontally and vertically along lines of constant NG 


and NS for each matrix coefficient. The Smoothed Dynamic 
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TABLE 1 


STEADY STATE CONVERGENCE MAP "A" MATRIX 
COEFFICIENT SIGN ERRORS 


GRID LOCATION MATRIX COEFFICIENT CORRECT SIGN 
(NG,NS) 
22000,1000 Al2 - 
22000,1500 Al2 - 
22000, 2000 Al2 - 
22000,2000 ro f 
22000, 2500 ne a 
23150,2695 iN, - 
23150,2695 eal + 
25000, 2500 Zl - 
25000, 3000 Al2 - 
25000, 3000 A21 + 
29100,500 A23 + 
29100,2950 Al2 - 
29100, 2950 A21 + 
3000,1000 A23 + 
32000,500 A23 + 
32000,1000 A23 i 


Transition Map of Figure 9 was formed by "eyeball smoothing" 
the data points from the Steady State Convergence Map. 
Trends were easily seen in the middle and right side of the 
Steady State Convergence Map, and these were used as models 


in the areas where sign discrepancies existed. The end 
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result was a grid of smoothed "A" matrices that would be 
the cornerstone of the dynamic simulation. 

A two-variable linear regression computer algorithm was 
used to obtain the coefficients of the best curve fit for 
each corresponding set of matrix coefficients. Table 2 is a 
summary of the six fits used. The A31, A32, and A33 entries 
of the "A" matrices were those that form the state variable 
corresponding to combustion energy and were always the same, 
hence no fit was required. The "B" matrix was constant at 
all values of NG and NS. 

The dynamic simulation of the plant was conducted using 
an IBM mainframe computer routine entitled Dynamic 
Simulation Language (DSL). The above described "A" matrix 
validation led to a method of successive linearizations to 
compute the values of NG and NS at time intervals of 0.001 
seconds during the dynamic trajectories being modelled. 

Strip chart data from actual plant runs was entered into 
the DSL code to provide the curve for model comparison. The 
initial and final plant setpoints were then specified, 
followed by the equations to compute the various matrix 
coefficients, the linearization equations, and the various 
output and graphing statements required. Figure 10 depicts 
Single and multiple linearizations plotted against data for 
a dynamometer speed versus time curve. Clearly, multiple 
linearizations were necessary to ensure the proper ending 


steady state values were reached. Appendix E 1s a copy of 


oo 


TABLE 2 


"A" MATRIX COEFFICIEND SCURRY Ena. 


MATRIX CURVE 
COEFFICIENT FIT EQUATION FORM 

All Exponential All = EXP(C1*NS+C2*NG+C3) 

Al2 Exponential Al2 = EXP(C1*NS+C2*NG+C3) 

A113 Exponential A1l3 = EXP(C1*NS+C2*NG+C3) 

Ron Exponential A21 = C1*NS“+C2*NS*NG 
+C3 *NG2+C4*NS 
+C5*NG+C6 

Bee Linear A22 = C1*NS+C2*NG+C3 

A23 Exponential A23 = EXP(C1*NS)+C2*NG+C3) 


the dynamic simulation program. Appendix F compares the "A" 
matrix coefficients of the Smoothed Dynamic Transition Map 
to those obtained by the dynamic simulation. 

Model validation was conducted in the region of the 
Smoothed Dynamic Transition Map known to be most reliable. 
Three runs were made across the map at constant NG speeds of 
23150, 29100, and 34900 rpm with varying NS speeds. The 
results are shown in Appendix G. All show excellent 
agreement between the model and data. A fourth validation 


was attempted vertically on the left side of the smoothed 


grid, starting at NG = 17000.0 rpm and NS = 500 rpm and 
ending at NG = 35500 rpm and NS = 748 rpm. The results 
obtained were less than satisfactory. However, this 


particular validation effort began and ended well outside of 


40 


0;00cE 


O°O00€ 


0°00Ge 


GO 0002 O'OOc) O'UQO0T 
(Ndu) GaadS LIVUS 


4] 


ONS 


N 


Al 


NEARI 


SINGLE LINEARIZATIO 


eee 


SLR TP Gil le Delis 


SE sas aa 





i 
O°00G 


Q°O 


35.0 


30.0 


IMD 
a 


20.0 


EC) 


10.0 


0.0 


(S 


TRI: 


EBrRreCcre OF Limearrzattrons om Ne 


Figune Or 


the region considered to be reliable on the "A" map, and 
progressed through the unreliable low dynamometer speed 
range. For these reasons the speed ranges chosen for the 
controller design and validation portion of this work were 
in the center of the Smoothed Dynamic Transition Map which 


was considered validated. 


42 


VI. CONTROLLER DESIGN AND COMPARISON 


The final segment of this work was concerned with 
controller design and comparison for overall performance and 
disturbance rejection qualities. Two controllers were 
designed, one a classical proportional integral controller 
(PI) and the other a modern linear quadratic regulator 
(LQR). Due largely to the fact that marine gas turbine 
controller designs are largely proprietary in nature, the 
design procedures used in this work had to be developed by 
the author. 

Deceleration was chosen as the design case because 
initial work clearly showed a much smaller stability margin 
than acceleration operations. Each controller was designed 
and validated for varying gas generator and shaft speed, but 
at constant dynamometer water volume. These specifications 
Simulate acceleration and deceleration modes at constant 
propeller pitch. Deceleration began at steady state speeds 
of NG = 34900 rpm and NS = 2000 rpm, and were subsequently 
slowed to NG = 23150 rpm and NS = 1500 rpm. The values were 
Simply reversed for the acceleration case. 

The control scheme shown in Figure 1 is currently 
employed by the U.S. Navy for control of General Electric 
IM-2500 gas turbine propulsion plants aboard DD-963 class 


destroyers. Note that the Navy uses a two loop approach to 
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control the gas turbine speed, propulsion shaft speed, and 
propeller pitch. Inputs to the plant were through an 
integrated throttle arrangement that schedules these three 
quantities for the desired ship speed. In order to 
demonstrate a similar control concept, the plant model used 
in the present work was chosen to closely emulate the 
general structure of Figure 1. Figure 11 is a diagram of 
this model and the control structure used for the PI design. 

The PI controller was designed first using a two step 
process. The inner loop consisted of a proportional 
control scheme to govern the gas generator. It was designed 
prior to the outer loop, using a cut and try process to 
select the appropriate gain (KPNG) to give a smooth non- 
oscillatory response. Steady state error was not a major 
consideration in the inner loop because the overall plant 
control was to be exerted on the propulsion shaft. Figures 
12 and 13 depict different responses for various choices of 
KPNG for deceleration and acceleration respectively. A gain 
of KPNG = 0.001 was chosen for the inner loop and held 
constant throughout the design of the outer loop. 

A proportional integral arrangement was employed for the 
outer NS control loop. This was chosen so that the steady 
state error associated with shaft response would be 
eliminated. Once again a cut and try procedure was used to 
decide the shaft proportional gain (KPNS) and integral 


(KINS) gain, with the criteria of smooth response driving 
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the selections. Figures 14 and 15 show plant response to 
various gains for the deceleration case in the gas generator 
and power turbine output shaft respectively. Gains of KPNS 
= 80.0 and KINS = 40.0 were selected for the outer control 
loop. 

The LQR controller was designed next using the control 
design software package MATRIXX. State space "A" and "B" 
matrices at the center of the smooth grid (NG = 25000 rpm, 
NS = 1500 rpm) were chosen to fix the design. In the LOR 
method, gains are sought to minimize a specified performance 
index "J" (or cost function) ff Ret. ela The performance 
index 1S expressed as an integral containing a function of 


the state variables and a function of the input variables: 


OO 


Cy 
Il 
or 


(eTQe + uTRu) dt (3) 


The "Q" and "R" matrices are symmetric weighting matrices 
that weight the states and inputs respectively. The 
designer chooses "Q" and "R", then computes the performance 
index which results in the LOR gains. Tradeoffs between "Q" 
and "R" weightings can be performed to achieve the best 
results. 

In the present work, the "Q" matrix was scaled so that 
each state matrix was making an approximately equal 
contribution to the response. The "R" matrix was chosen by 


a cut and try process, and was designed to minimize the 
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plant fuel input. The matrices were adjusted until the most 
acceptable response was obtained. This particular LOR de- 
sign is strictly a proportional regulator and has no inte- 
gral action to remove steady state error. As a result there 
is some steady state error in both NG and NS. It was possi- 
ble to eliminate or nearly eliminate this error in either NG 
or NS, but at the expense of the other variable. The chosen 
design exhibits small steady state errors in both NG and NS. 

Figures 16, 17, 18, and 19 are comparisons of the 
deceleration and acceleration validations of the PI and LOR 
controllers for both the gas generator and power turbine 
shafts. The PI controller provides a smoother response in 
both NG and NS deceleration curves, while LOR reaches steady 
state in less time. The acceleration validation shows the 
LOR controller providing a smoother, quicker response in NS 
and a quicker response in NG. In actuality, all responses 
are comparable for both controllers. 

Disturbance rejection for both controllers was analyzed 
by subjecting each one to a cyclic torque disturbance 
Simulating sea state oscillations. A sine function load was 


used that provided a ten second period: 


Qr = 20.0 sin(t™/5.0) (4) 


The controllers were set to maintain steady gas generator 


and shaft speeds of NG = 17642 rpm and NS = 1500 rpm 
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respectively, and the Q;, function applied through the torque 
input to the models. Each model was run for 35 seconds and 
the data recorded. Peak-to-peak values of NG and NS were 
obtained, as well as fuel consumed by the plant for both 
control designs. Table 3 compares the results for two LOR 
designs to the PI design and graphically illustrates the 
effects that the before mentioned tradeoffs can have on the 
LQR design. Of particular interest are the NG peak-to- peak 
values and the fuel consumed values. The smaller NG peak- 
to-peak values indicate better disturbance rejection for the 
NG output, which translates to less wear on critical gas 
turbine components such as variable stator vanes. This in 
turn has the potential to decrease maintenance downtime as 
well as mean time between failures of complex mechanical 
components. Although the difference in fuel consumed may be 
small, the potential exists for substantial fleetwide 
reductions in fuel costs when this figure is applied to the 


amount of fuel burned annually by all gas turbine ships. 
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TABLE 3 


CONTROLLER DESIGN COMPARISON 


nct nst FUEL2 
LOR #1 3320.0 DOB oN Ong 4136 
LOR #2 3688.0 125.8 0.74208 
PI 4233.0 74.5 Oar 4ea7 


lpeak-to-Peak Values, rpm 


Total lb, consumed in 35 seconds 
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VII. CONCLUSIONS AND RECOMMENDATIONS 


Two non-proprietary control design methods have been 
developed for marine gas turbine propulsion systems, one a 
Classical PI controller, the other a modern LOR controller. 

Comparable performance has been demonstrated between the 
PI and LOR controllers. 

A new gas turbine simulation has been developed that 
allows real time or near real time computation of performn- 
ance. This simulation has immediate applications in model 
based control and plant health monitoring. 

To increase confidence in the plant model, the Steady 
State Convergence Map should be reconstructed from data 
obtained from actual plant runs at designated points 
throughout the operating envelope. This would eliminate the 
problems of multiple root convergence in the steady state 
computer simulations, and provide a more accurate data base 
for the curve fits required by the dynamic simulation. 

Proportional LQR design should be further developed to 
account for important non linearities in the gas generator 
and controllable reversible pitch propeller, as well as 
limiting/alarm conditions that exist in fleet marine gas 
turbine propulsion plants. 

An integral LOR controller should be investigated. This 


type of controller has the potential to offer better 
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performance tradeoffs in the areas of gas generator speed 
control, shaft speed control, and in the reduction of fuel 
consumption. 

Once the LOR controller is perfected, the idea of an LOR 
Gain Map similar to the Smoothed Dynamic Transition Map 
should be investigated. This concept has the potential to 
Maximize the benefits of the LQR controller by computing the 
best gain values at each time step as the plant progresses 
through a particular dynamic trajectory. 

This technology should be implemented at the Naval 
Postgraduate School to quantify real improvements’) and 


justify larger scale development. 
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ROUTINE VRD 


By 
V.A. STAMMETTI 


DL. ssi 


THIS ROUTINE PROVIDES THE FIRST OF A THREE STEP 
PROCESS TO OBTAIN STEADY STATE CONVERGENCE AT A POINT 
IN THE OPERATING ENVELOPE OF THE BOEING 501-6A GAS 
TURBINE INSTALLED AT THE NAVAL POSTGRADUATE SCHOOL. 
ROUTINE VRD IDENTIFIES THE SEARCH RANGES FOR THE 


VARIABLES MF, P2, AND P4. 
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seveaedesesedesedevevededesededededevedevevetedsdedcsesvevededvevedssededesvesedcdetevedededescsevevevetesvevedede deve devedetesesevedve 


REAL NG,NS,MF,MA,MAF ,MFG,MFI,MFHIGH ,MFLOW ,MFSAVG 


INPUT THE INITIAL GAS GENERATOR SPEED AND DYNO SPEED. 


WRITE( 6,2) 
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FORMAT(/,3X,' INPUT INITIAL GAS GENERATOR SPEED,''NG”. 


READ(5,*) NG 
WRITE(6,*) NG 


WRITE(6,3) 


FORMAT(/,3X,' INPUT INITIAL DYNO SPEED,"NS”. ') 


READ(5,*) NS 


WRITE(C6,*) NS 


WRITE(9,*) ' NG=' ,NG 


WRITE(9,*) ' NS=',NS 


WRITE(6,4) NMF 
FORMAT( /,3X,' INPUT NUMBER OF FUEL GUESSES,'NMF". ' ) 
READ(5,*) NMF 


WRITE(6,*) NMF 


WRITE(6,5) NP2 

FORMAT(/,3X,' INPUT NUMBER OF P2 GUESSES,''NP2”. ' ) 
READ(5,*) NP2 

WRITE(6,*) NP2 
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WRITE(6,6) NP4 
FORMAT(/,3X,' INPUT NUMBER OF P4 GUESSES,"NP4".') 
READ(5,*) NP4 


WRITE(6,*) NP& 


CALCULATION OF INITIAL FUEL GUESS 


CALL NGNSMFCNG,NS ,MFT) 


VARIABLE INITIALIZATION 


MEG. = Mh e200 


RNMF = NMF 


NP2 


RNP2 
RNP4 = NP4 


DMF 


40.0 / RNMF 


DP2 25.0 / RNP2 


DP4 = 5.0 / RNP4 


QCONV = 10.0 
P2CONV = 0.5 
P4CONV = 0.5 


MFLOW = MFI 

MFHIGH = MFG 
P4LOW = 100.0 
P4HIGH = 1.00 
P2LOW = 100.0 
P2HIGH = 1.00 


ATEST = 1000.0 
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P2SAVC = 0.0 
P2SAVG = 0.0 
P4SAVC = 0.0 
P4SAVG = 0.0 
MFSAVG = 0.0 

MF LOOP 

DO 100 J = 1,NMF 


MFG = MFG + DMF 
P2SAVE = 0.0 
P2G = 20.0 
WRITE(6,*) 


WRITE(6,*) '‘MFG=' ,MFG 


P2 LOOP 


DO 200 K = 1,NP2 

P2G = P2G + DP2 

CALL SUBQC(NG,P2G,QC) 
CALL SUBT2(NG,P2G,T2) 
CALL SUBMA(NG, P2G,MA) 
QSAVE = 0.0 

P4SAVE = 0.0 


P4G = 20.0 


P4& LOOP 
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10 


[> 


20 


DO 300 L = 1,NP4 


P4G = P4G - DP4 


TORQUE CONVERGENCE 


CALL SUBQHT(NG ,MA,T2,MFG,P4G,QHPT) 


DELQ = QC - QHPT 


QTEST = DELQ * QSAVE 


QSAVE 


DELQ 


IF @@REST. GE. 0.0). GO TO 300 


P4 CONVERGENCE 


MAF = MA + MF 

CALL SUBT4(NG,MA,T2,MFG, P4G,T4) 
CALL SUBP4(MAF ,T4,NS, P4C) 

DELP4 = P4C - P4G 

P4TEST = DELP4**P4SAVE 

P4SAVE = DELP4 

IF(P4TEST. GE.0.0) GO TO 300 


IFC ABS(DELP4). GT. P4CONV) GO TO 300 


P2 CONVERGENCE 


CALL SUBP2(NG,MA,T2,MFG,P4G,P2C) 
DEER? = P20 = -P2G 
P2TEST = DELP2*P2SAVE 


P2SAVE 


DELP2 


30 


35 


Tet? 20ES1. GE. 0.0) GO 10 300 


IF(MFG. LT. MFLOW) MFLOW = MFG 
TFCMEG.GI.MFHIGH) MFHIGH = MFG 
IF(P4G. LT. P4LOW) P4LOW = P4G 
IF( P4G. GT. P4HIGH) P4HIGH = P4G 
IF(P2G. LT. P2LOW) P2LOW = P2G 
IF( P26. GT. P2HIGH) P2HIGH = P2G 


DELSUM = (DELQ**2) + (DELP2**2) + (DELP4**2) 


IF(DELSUM. GE. ATEST) GO TO 300 


P2SAVC = P2C 
P2SAVG = P2G 
P4SAVC = P4C 
P4SAVG = P4G 
MFSAVG = MFG 
DELQG = DELQ 


DELP2G = DELP2 
DELP4G = DELP4 


ATEST = DELSUM 


WRITE(6,*) ' CONVERGENCE OBTAINED AT' 


WRITE(6,*) ‘ P2C=',P2SAVC 
WRITE(6,*) ' P2G=' ,P2SAVG 
WRITE(6,*) ' P4C=',P4&SAVC 
WRITE(6,*) ' P4G=',P4SAVG 
WRITE(6,*) ' MF=' ,MFSAVG 

WRITE(6,*) ‘ DELQ=' ,DELQG 
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300 


200 


100 


WRITE(6,*) ' DELP2=' ,DELP2G 


WRITE(6,*) ' DELP4=' ,DELP4G 


CONTINUE 


CONTINUE 


CONTINUE 


WRITE(6,*) ' CONVERGENCE OBTAINED AT' 


WRITE(6,*) ‘ P2C=' ,P2SAVC 
WRITE(6,*) ‘ P2G=',P2SAVG 
WRITE(6,*) ' P4C=' ,P4SAVC 
WRITE(6,*) ' P4G=' ,P4SAVG 
WRITE(6,*) ' MF=',MFSAVG 
WRITE(6,*) ' DELQ=' ,DELQG 
WRITE(6,*) ' DELP2=' ,DELP2G 
WRITE(6,*) ' DELP4=' ,DELP4G 


WRITE(6,*) ' MFLOW=' ,MFLOW 
WRITE(6,*) ' MFHIGH=' ,MFHIGH 
WRITE(6,*) ' P4&LOW=' ,P4LOW 
WRITE(6,*) ' P4HIGH=' ,P4HIGH 
WRITE(6,*) ' P2LOW=' ,P2LOW 


WRITE(6,*) ' P2HIGH=' ,P2HIGH 
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900 


WRITE(9,*) 
WRITE(9 ,*) 
WRITE(9 ,*) 
Whepec9 ,*) 
WRITE(9 ,*) 
WRITE(9,*) 
Wier CS), 7 ) 
WE 9,7) 
WRITECS,*) 


WRITE(9,**) 
WRITE(9,*) 
WRITE(9,**) 
WRITE(9,**) 
WRITE(9,**) 


WRITE(9,*) 


STOP 


END 


f 


t 


1 


CONVERGENCE OBTAINED AT' 
P2C=' ,P2SAVC 
P2G=' , P2SAVG 
P4C=' , P4SAVC 
P4G=' , P4SAVG 
MF=' ,MFSAVG 
DELQ=' , DELQG 
DELP2=' , DELP2G 


DELP4=' , DELP4G 


MFLOW=' , MFLOW 
MFHIGH=' ,MFHIGH 
P4LOW=' , P4LOW 
P4HIGH="' , P4HIGH 
P2LOW=' , P2LOW 


P2HIGH=' , P2HIGH 
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APPENDIAX B 


Vee eres Tedd eT Tees Be PRET E TE TE TEDE BETES TEES TEETER ERGOT TC DTEREERERERERELIEREREITC 
ROUTINE SC ze 

¥ 

BY v 

V.A. STAMMETTI a 

De Leer tL % 

THIS ROUTINE PROVIDES THE SECOND OF A THREE STEP % 
PROCESS TO OBTAIN STEADY STATE CONVERGENCE AT A POINT * 
IN THE OPERATING ENVELOPE OF THE BOEING 3502-64 GAS Le 
TURBINE INSTALLED AT THE NAVAL POSTGRADUATE SCHOOL. % 
ROUTINE SC CONVERGES THE VARIABLES MF, P2, AND P4. co 
VETTE TET S TOTE TOTES OTE TEETER ETRE EB RRBRERE ERE SREREREBREEREERREREERRREEREE 


seveves's 


C 
c 


REAL*8 NG,NS,MF,MA,MAF ,MFG,MFI ,MFSAVG,MFHIGH ,MFLOW ,MFL,RNMF , 
RNP2 , QCONV , P2CONV,, P4CONV, P2LOW , P2HIGH, P4J.OW, P4HIGH, ATEST, 
P2SAVE , P2SAVC, P2SAVG , P4SAVE , P4SAVC, P4SAVU , UMF , DP2 ,DP4, RNNMF, 
DMF2,P2L,RNNP2 ,DP22,P4H,RNNP4 ,DP42,P2G,QC,T2,QSAVE, P4G, PAC, 
P2C ,QHPT,DFLQ,QTEST,T4, DELP4, P4TEST, P2TEST, DELP2 ,DELSUM, 


RNP4 , DELQG , DELP4G , DELP2G 


INPUT THE INITIAL GAS GENERATOR SPEED AND DYNO SPEED. 


WRITE(6, 2) 


FORMAT(/,3X,'INPUT INITIAL GAS GENERATOR SPEED,''NG". ') 


READ(5,*) NG 


WRITE(6,*) NG 
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WRITE(6,3) 


FORMAT(/,3X,' INPUT INITIAL DYNO SPEED,'"NS". ') 


READ(5,*) NS 
WRITE(6,*) NS 


WRITE(9,*) ' NG=',NG 


WRITE(9,*) ' NS=' ,NS 


WRITE(6,4) NMF 
FORMAT(/,3X,'INPUT NUMBER OF FUEL GUESSES,"NMF". ') 
READ(5,*) NMF 


WRITE(6,*) NMF 


WRITE(6,5) NP2 
FORMAT(/,3X,' INPUT NUMBER OF P2 GUESSES,"NP2". ') 
READ(5,*) NP2 


WRITE(6,*) NP2 


WRITE(6,6) NP4& 
FORMAT( /,3X,' INPUT NUMBER OF P4& GUESSES,"NP4". ') 
READ(5,*) NP4 


WRITE(6,*) NP4 
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WRITE(6,7) NRNO 
FORMAT(/,3X,' INPUT THE NUMBER OF THIS REFINEMENT,"NRNO". ') 
READ(5,*) NRNO 
WRITE(6,%*) NRNO 


VARIABLE INITIALIZATION 


RNMF = NMF 
RNEZ = NEZ 
RNP4 = NP4 


QCONV = 10.0 


P2CONV Cx 
P4CONV = 0.5 

MFLOW = 86. 0000000 
MFHIGH = 90. 00000 
P4LOW = 15. 000000 
P4HIGH = 16. 0000000 
P2LOW = 22. 0000000 
P2HIGH = 23.000000 
ATEST = 1000.0 
P2SAVC = 0.0 
P2SAVG =0.0 

P4SAVC =0.0 


P4SAVG 


lI 
© 
© 


MFSAVG = 0.0 


DMF = 40.0 / RNMF 
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DP2 = 25.0 / RNP2 
DP4 = 5.0 / RNP4 
MFL = MFLOW - DMF 


WRITE(6,*) 'MFL=' ,MFL 

NMF = (MFHIGH - MFL) / DMF 
NNMF = NMF * 5 

RNNMF = NNMF 

DMF2 = DMF / 5.00 

NNMF = NNMF + 10 
WRITE(6,*) 'NNMF=' ,NNMF 


WRITE(6,*) 'DMF2=' ,DMF2 


P2L = P2LOW - DP2 
WRITE(6,*) 'P2L=' ,P2L 

NP2 = (P2HIGH - P2L) / DP2 
NNP2 = NP2 * 10 

RNNP2 = NNP2 

DP22 = DP2 / 10.0 

NNP2 = NNP2 + 10 
WRITE(6,*) 'NNP2=' ,NNP2 
WRITE(6,*) 'DP22=' ,DP22 
P4H = P4HIGH + DP4 

NP4 = (P4H - P4LOW) / DP4 
NNP4 = NP4 * 10 

RNNP4 = NNP4 

DP42 = DP4 / 10.0 

NNP4 = NNP4 + 10 


WRITE(6,*) 'NNP4=' ,NNP4 


al 


C2. (C2) 3G 


WRITE(6,*) ‘DP42=' ,DP42 


REINITIALIZE VARIABLES 


DMF = (MFHIGH - MFLOW)/RNMF 
DP2 = (P2HIGH - P2LOW) / RNP2 
DP4 = (P4HIGH - P4LOW) / RNP4 


MFHIGH = 20.0 


MFLOW = 500.0 


P2LOW 800), 0 
P2HIGH = 1.0 
P4LOW = 300.0 


P4HIGH = 1.0 


MFG<= EL ="piire 
MF LOOP 


MFG = MFLOW - DMF 


DO 100 J i NNIE 
DO 100 J = 1,NMF 
MFG = MFG + DMF2 
MEG — MEG +R: 
P2SAVE = 0.0 
P2G°="P2b = DPZ2 
P26 — =P 2LOW F-eDer 


WRITE(6,*) 


oT Or Pee 


WRITE(6,*) 'MFG=' ,MFG 


PZ SE@OOP 


DO 200 K = 1,NNP2 


DO 200 K = 1,NP2 
Pz2G — P2G 4 >0R2 
E2G = P26 Drz2 


CALL SUBQC(NG,P2G,QC) 
CALL SUBT2(NG,P2G,T2) 
CALL SUBMA(NG,P2G,NA) 
QSAVE = 0.0 

P4SAVE = 0.0 

P4G = P4H + DP42 

P4G = P4HIGH - DP4 


WRITE(6,*) 'P4G=' ,P4G 


P4& LOOP 


DO 300 L = 1,NNP4 


P4G = P4G - DP42 


DO 300 L = 1,NP4 


P4G = P4G + DP4 


TORQUE CONVERGENCE 


CALL SUBQHT(NG,MA,T2,MFG,P4G,QHPT) 


DEE@r= QC = QHPT 
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10 


14 


i 


20 


DELQ * QSAVE 


QTEST 


QSAVE = DELQ 

WRITE(6,*) '‘QTEST=' ,QTEST 
IF(QTEST.id. 0. 0)GOe10 10 
GO TO 300 


IFC ABS( DELQ). GT. QCONV) GO TO 300 


P4 CONVERGENCE 


MAF = MA + MF 

CALL SUBT4(NG,MA,T2,MFG,P4G,T4) 
CALL SUBP4(MAF,T4,NS,P4C) 

DELP4 = P4C - P4G 

P4TEST = DELP4*P4SAVE 

P4SAVE = DELP4 

WRITE(6,%*) 'P4TEST=' ,P4TEST 
IF(P4TEST. LE. 0.0) GO TO 20 

GO TO 300 


IFC ABS( DELP4).GT. P4CONV) GO TO 300 


P2 CONVERGENCE 


CALL SUBP2(NG,MA,T2,MFG, P4G, P2C) 


DEEP2 = P20 = 3226 


P2TEot = DELP2*P2SAVE 


P2SAVE DEBEZ 
WRITE(6,*) ‘P2TEST=' ,P2TEST 


TIPCP2ZIEST. LE. Os 0) GG eee 
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24 
25 


30 


aD 


GO TO 300 

TE GAZS@DELP2). GT. P2CONV) GO T@_300 
IF(MFG. LT. MFLOW) MFLOW = MFG 
IF(MFG. GT. MFHIGH) MFHIGH = MFG 
IF(P4G. LT. P4LOW) P4LOW = P4G 
IF(P4G. GT. P4HIGH) P4HIGH = P4G 
IF(P2G. LT. P2LOW) P2LOW = P2G 
LhQe2G. Gl. P2iGH) ePZuIGi =" P26 


WRITE(6,*) ‘CONVERGENCE HERE’ 


DELSUM = (DELQ**2) + (DELP2**2) + (DELP4**2) 


PEGDELSUM.GE. ATEST)mGOsTOs300 


P2SAVC Ee 
P2SAVG = P2G 
P4SAVC = P4C 
P4SAVG = P4G 
MFSAVG = MFG 
DELQG = DELQ 
DELP2G = DELP2 
DELP4G = DELP4 


ATEST = DELSUM 


WRITE(6,*) ' CONVERGENCE OBTAINED AT' 


WRITE(6,*) ' P2C=',P2SAVC 
WRITE(6,*) ' P2G=',P2SAVG 
WRITE(6,*) ' P4C=',P4SAVC 
WRITE(6,*) ' P4G=',P4SAVG 


7S 


300 


WRITE(6,*) 
WRITE(6,*) 
WRITECG,*) 
WRITE(6,*) 
WRITE(6,*) 


WRIDE CR) 
WRITEC9. 
WRITE CoE) 
WRITEC9,*) 
WRITEC9,*) 
WRITE(9,*) 
WRITECOS) 
WRITE(C9,*) 


WRITE(9,*) 


WRITE(9,**) 
WRITE(9,*) 
WRITE(9,**) 
WRITE(9,**) 
WRITE(9 ,**) 
WRITE(9 ,**) 
WRITE(9,**) 
WRITE(9,*) 
WRITE(9,**) 


— 


t 


MF=' ,MFSAVG 
DELQ=' , DELQG 
DELP2=' , DELP2G 
DELP4=' , DELP4G 


NUMBER OF REFINEMENTS=' ,NRNO 
CONVERGENCE OBTAINED AT' 
P26— 2 P25A06 

P2G=' ,P2SAVG 

P4C=' , P4SAVC 

P4G=' , P4SAVG 

MF=' ,MFSAVG 

DELQ=' , DELQG 

DELP2=' , DELP2G 


DELP4=' , DELP4G 


MFLOW=' , MFLOW 
MFHIGH=' ,MFHIGH 
P4LOW=' , P4LOW 
P4HIGH=' , P4HIGH 
P2LOW=' , P2LOW 


P2HIGH=' , P2HIGH 


P4 = P&G + DP42 


CALL PART(CA,B,MA,P2G,T2,MFG,NG,NS,P4G,T4, MAF ,WW) 


CONTINUE 
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Gy tay ta OD 


200 


100 


900 


CONTINUE 


CONTINUE 


WRITE(6 ,*) 
WRITE(6,*) 
WRITE(6,*) 
WRITE(6,*) 
WRITE(6,**) 
WRITE(6,*) 
WRITE(6,*) 
WRITE(6,**) 
WRITE(6 ,*) 


WRITE(6,*) 
WRITE(6,*) 
WRITE(6,*) 
WRITE(6,*) 
WRITE(6,*) 
WoldECo,*) 


STOP 
END 


CONVERGENCE OBTAINED AT' 
P2C=' ,P2SAVC 

P2G=' ,P2SAVG 
P4C=' , P4SAVC 
P4G=' , P4SAVG 

MF=' ,MFSAVG 
DELQ=' , DELQG 
DELP2=' , DELP2G 


DELP4=' , DELP4G 


MFLOW=' ,MFLOW 
MFHIGH=' ,MFHIGH 
P4LOW=' , P4LOW 
P4HIGH=' , PGHIGH 
P2LOW=' ,P2LOW 
P2HIGH=' , P2HIGH 


a 


APPENDIX C 


tevededtevedesededededtedesevededeteckdtevetetkdestelededesedesiededeketeteiehhetecdehiteksesesedtetedidtetcdekhhkheeReikKsk 


= Hs 
* ROUTINE SSVPD * 
: ve 
* MODIFIED BY ae 
V.A. STAMMETTI ¥ 

* FROM A SUBROUTINE BY ve 
* V. J. HERDA or 
LL x 
we THIS ROUTINE IS THE THIRD AND FINAL OF A THREE STEP % 
* PROCESS TO OBTAIN STEADY STATE CONVERGENCE AT A POINT * 
IN THE OPERATING ENVELOPE OF THE BOEING 501-6A GAS we 

* TURBINE ENGINE INSTALLED AT THE NAVAL POSTGRADUATE * 
* SCHOOL. ROUTINE SSVPD PROVIDES THE STEADY STATE * 
* VALUES FOR ALL NECESSARY PLANT VARIABLES, AS WELL AS * 
* THE STATE SPACE "A" AND "B" MATRICES. ce 
ry se 


teeetotkdeleletcticlhkiohkhkhewcetokvtkionhethkkidielhkkhnhkKhrThkkeakReikakekReRKRR RR: 


REAL*8 NG,NS,MF,MA,MAF,NSO,NGO,MFO,MAFO,MAO,MFDEL,MFU,MFL, 


1 MIR,MF2,MFMIN,MERR,MAC,AFC 72,14) 24¢)P2G5 am 
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LROOT , HROOT , P4G1 , P4G2 , P4TEST,AA,BB,CC,AA1,BB1,CCl1, 
SCALE , P2DEL, P4CONV, FO 
DIMENSION AC3,3),BC3,2) 


INPUT THE INITIAL GAS GENERATOR SPEED AND DYNO SPEED. 


WRITE(6, 2) 
FORMAT(/,3X,' INPUT INITIAL GAS GENERATOR SPEED,'"'NG". ' ) 


READ(5,*) NG 


WRITE(6,*) NG 


WRITE(6,3) 


FORMAT(/,3X,' INPUT INITIAL DYNO SPEED,"NS”.') 


READ(5,*) NS 


WRITE(6,*) NS 


WRITE(6,4) P2 
FORMAT(/,3X,' INPUT PRESSURE,'P2". ') 
READ(5,*) P2 


WRITE(6,*) P2 


WRITE(6,5) NP4& 


FORMAT(/,3X,' INPUT PRESSURE,''P4". ') 
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READ(5,*) P4 


WRITE(6,*) P4 


WRITE(6,6) MF 
FORMAT(/,3X,' INPUT FUEL,''MF”. ') 
READ(5,*) MF 


WRITE(6,*) MF 


PARAMETER CALCULATIONS 


CALL SUBMA(NG,P2,MA) 

CALL SUBT2(NG, P2,T2) 

CALL SUBT4(NG,MA,T2,MF,P4,T4) 
CALL SUBQC(NG, P2,QC) 

CALL SUBQHT(NG,MA,T2,MF,P4,QHPT) 
MAF = MA + MP 


CALL SUBQFT(MAF ,T4,NS,QFPT) 


QD = QFPT 
CS = 1. 1920e-5 
C3 =-4,0E-6 


C& = -20.0 + C3*NS*NS 


MIR = (QD - C4) / (C5*NS*NS) 


IF(MIR: LT. 0. 0) 9MIR ="0N0 


WW = MIR**( 1.0/1. 3) 


COMPUTATION OF THE STATE SPACE MATRIX COEFFICIENTS 


CALL PART(A,B,MA,P2,T2,MF,NG,NS,P4,T4,MAF,WW) 
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Co Ga Ga 


WRITE(6,*) 
WRITE(6,*) 
WRITE(6,*) 
WRITE(6,*) 
WRITE(6,*) 
WRITE(6,*) 
WRITE(6,**) 
WRITE(6,*) 
WRITE(6,”*) 
WRITE(6,*) 
WRITE(6,*) 
WRITE(6,*) 
WRITE(6,*) 


WRITE(6,*) 


WRITE(9,*) 
WRITE(9,*) 
WRITE(9,*) 
WRITE(9,*) 
WRITE(9,*) 
WRITE(9,*) 
WRITE(9,*) 
WRITE(9,*) 


WRITE(9,*) 


NG = ',NG 
NS = ',NS 
MF = ' MF 
P2 = ",P2 
P4 = ',P4 
2 2 
T4 = ',T4 
MA = * ,MA 
MAF = ' ,MAF 
ac = ',Qc 
QHPT = ' ,QHPT 
QFPT = ' ,QFPT 
QD = ',QD 
WW = ',WW 
NG = ',NG 
NS = ',NS 
MF = ' MF 
Rea Re 
P4 = ',P4 
T2 = "72 
T4 = ',T4 
MA = "= MA 
MAF = ' ,MAF 
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WRITE(9,*) 
WRITE(9,*) 
WRITE(9,*) 
WRITE(9,*) 


WRETEC ee) 


900 STOP 
END 


C 


ac = *,Qc 
QHPT 


' ,QHPT 
QFPT = ',QFPT 
QD = ',QD 


Cre tererere vend devevevedede te deve ve de Ts te ve ve ess ve ve ys Ve Ne Te Nee Ve Me Ve Te Te Fe Fe Ve Ie Ve He Fe Ne Me Me Ve Me Me Vee Ve Me TeV Te Ne Me 


SUBROUTINE PART(A,B,MA,P2,T2,MF,NG,NS,P4,T4,MAF ,WW) 


Crevededsaevenevevede de ue ve vee dene sede ve seve ve deve Te Fe se 1 Tete Be Me Fe Fe Me Me Te Me Mee Tes es Te TENE Tee Me Tee Ie Tee TEN Te 


C 


C THIS SUBROUTINE CALCULATES THE ELEMENTS OF THE 'A' AND 'B' MATRICES 


C IN THE STATE SPACE EQUATION: 


c 
E XDOT = A*X + BYU . 
C 
c 
C COMMON QC,NG,P2,QH,MA,T2,MF,P4,QF ,MAF,T4,NS,QD,WW 
DIMENSION A(3,3),B(3,2) 
REAL NG,NS,MF,MA,MAF,JG,JD,DENOM1,DENOM2,DENOM3 
C 
C 
C 
JG = 0.009525 * 2 * 3.14159 / 60.0 
JD = 0.6738 * 2 * 3.14159 / 60.0 


C 


C CALL SUBROUTINES TO GET PARTIAL DERIVATIVES. 


C 
CALL DMA(NG, P2,DMADNG, DMADP2) 
CALL DT2(NG,P2,DT2DNG,DT2DP2) 
CALL DQC(NG, P2,DQCDNG , DQCDP2) 
CALL DP2(NG,MA,T2,MF,P4,DP2DNG , DP2DMF , DP2DMA, DP2DT2 , DP2DP4) 
CALL DT4(NG ,MA,T2,MF,P4,DT4DNG , DT4DMF , DT4DMA , DT4DT2 , DT4DP4) 
CALL DQHT(NG,MA,T2,MF,P4,DQHDNG, DQHDMF , DQHDMA , DQHDT2 , DQHDP4 ) 
CALL DP4(MAF ,T4,NS,DP4DNS , DP4MAF , DP4DT4) 
CALL DQFT(MAF,T4,NS,DQFDNS , DQFMAF , DQFDT4) 
CALL DQD(NS, WW, DQDDNS , DQDDWW) 

C 


C COMPUTE THE COEFFICIENTS OF THE STATE SPACE EQUATIONS (I.E. THE 


C ELEMENTS OF THE 'A' AND 'B' MATRICES). 


C 
J1 = DQHDMA 
J2 = DQHDMF 
J3 = DQHDT2 
J4 = DQHDP4 
JS = DQHDNG 
El = DQCDP2 
E2 = DQCDNG 
C1 = DMADP2 
C2 = DMADNG 
Dl = DT2DP2 
D2 = DT2DNG 
Al = DP4MAF 
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(H5 + H1"C2 + H3*D2) / DENOM1 


A2 = DP4DT4 

A3 = DP4DNS 

Hl = DP2DMA 

H2 = DP2DMF 

H37= DR2ZDT2 

H4 = DP2DP4 

H5 = DP2DNG 

G1 = DT4DMA 

G2 = DT4DMF 

G3 = DT4DT2 

G4 = DT4DP4 

G5 = DT4DNG 

Bl = DQFMAF 

B2 = DQFDT4 

B3 = DQFDNS 

Zl = B2*G3*)2. 4 8265 
Z2 = Bl + B2*G2 
Z3 = Bl + B2*G1 
Z4 = B2*G3*D1 

Z5 = B2*G4 

20 = Globe; 3562 
Zi = 2&4 + Z3*C1 
DENOM] = 1.0 = HICSS 
Yl = 

Y2 = H2 / DENOM1 
Y3 = H4 / DENOMI1 


DENOM2 = 1.0 - A2*G4 


Y4 = (A2*G5 + A1*C2 + A2*G3*D2 + A2*G1*C2) / DENOM2 
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= (A1*C1 + A2*G1*C1 + A2*G3*D1) / DENOM2 


¥5.=9A305/ DENOM2 

Y6 = (Al + A2*G2) / DENOM2 
nae 

DEWOMo = 1e0= Y7*Y5 

Y8 = (Y4 + Y7*Y1) / DENOM3 
Y9 = Y5 / DENOM3 


Y10 = (Y6 + Y7*Y2) / DENOM3 


Z8 

29 

Z10 
nell 
eZ 
was: 
Y14 


n15 


Z12 


Z13 


Z6 + Z7*Y1 + Z5*Y8 + Z7*Y3*Y8 
fore 2) *X 5°19 + BS 


Gow CN 2 erie 1 Oma 7 * YY 3°°Y 10 


—ageor= E2 + giee2-+ J3eD)2 


Seo Clete JSD l= El 


=e elite Yl 228 ¥ 1 
= a2 + Y12*Y2 
= J& + Y12*Y3 
= Y13 + Y15*Y8 
= Y15*Y9 


Sey i4e YI5“Y10 


FINAL FORM OF THE 'A' AND 'B' MATRICES. 


f NOTE ! 


ELEMENTS A33 AND B31 ARE NOT COMPUTED HERE BUT WERE 


DETERMINED EXPERIMENTALLY FROM GAS TURBINE TEST DATA. 


FOR ACCELERATIONS USE: 


A33 


B31 
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FOR DECELERATIONS USE: 


A33 = -0. 87 
Bole s0557 

All = ZakiyoG 

Al2 = Zi iG 

A13 = 213/JG 

A211 = 76730 

A22 = 29/JD 

A23 = 210/JD 

A31 = 0.0 

A32 ="0,0 

A33 = -10.0 

Bll = 0.0 

Bi22— "00 

B21 = 0.0 

B22 =) 10D 

B319 =. 1000 

B32 = 0.0 

WRITE(9,*) '"A" AND "B' MATRICES FOR NG = ',NG 

WRITE(9,*) ' AND NS = ',NS 

WRITE(9,*) ' ' 

WRITE(9,*) ' ' 

WRITE(9,%*) | All = A 

WRITE(9,*) ' AOD =o ee 

WRITE(9,*) ' Al3 = ‘gage 
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WRppEC9 ,* ) 
WRITE(9,*) 
WRITE(9,*) 
WRUGEC IS) 
WRITEC 9. *) 
WRITE(9,*) 
WRITE CIS) 
WRITE(9 ,*) 
WRITEC9,*) 
WEDTE CG: 2°) 
Whoo, *) 
WRITEC 9,7) 


RETURN 
END 


A21 
A22 
A23 
A31 
A32 
A33 


pad 
Ba 
B22 
B31 


BeZ 


A21 
A22 
A23 
A31 
A32 
A33 
Bll 
Ba 
B21 
B22 
Bou 
BZ 
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APPENDIX D 


Sevededetedededesede dele KK RT CTE eT eC To TF o TEV FCS SSeS TTS EM TETE TA TERT ee te sede sere ae 
oe ve 
ve THE FOLLOWING SUBROUTINES WERE WRITTEN BY V.J. HERDA ¥ 
* AND ARE USED IN VARIOUS COMBINATIONS BY THE ROUTINES VRD, * 
* §C, AND SSVPD. ¥ 


ate 
ee 
seseseseslcvecedeveseveveveslevevosedesetededevesesedeteteseteseestedesetetetekedeevekscsesksesedesestickdekacsetesese 


Seacseseseskoe kak skoksesdeokakakeseskakscseveseakakeaese acdsee skdeveveskcsdcakseskvesede akeakeseskeake sede ese ake ake ake okakvsevese 


SUBROUTINE NGNSMF(X1,X2,BR) 


Cavdedededodedevede de deat ae voce vee te oe te de Fe Fe Fe Fe Fe eM ETE Fe Te Tee Fee Te TE Ie Te IE TE TE TE Te IE He FEE Te IE HEE HIE Te AE IE IE 


C THIS SUBROUTINE PRODUCES AN INITIAL "GOOD GUESS" FOR 'MF' 


C BASED ON THE SPECIFIED INPUTS, 'NG' AND 'NS'. 


DIMENSION XS )5CC21)5 20>) Xk) 


XR( 1) X1 


XR(2) = X2 


C COEFFICIENTS OF THE QUADRATIC CURVE FIT. 


C( 1)= 1.982237 


G@e)— 02461500 


8§ 


C 


C 


C 


C( 3)= -5. 147902E-02 
C( 4)= -1. 884269 
C( 5)= -9.572456E-02 


Coe )= 0. 7639415 


SCALING FACTORS. 
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Z(1) =  36000.0 
Z(2) = 3000.0 
Z(3) = 240.0 
NIND = 2 


DO 686 I = 1,NIND 
I yA 


CONTINUE 


CONSTRUCT THE COMPLETE QUADRATIC EQUATION. 


71 


70 


B 


0 
K=0 
DO 70 J = 1,NIND 
DO 71 I = J,NIND 

K = K+1 
B = BtC(K)*X(J)*X(I) 
CONTINUE 


CONTINUE 


89 


72 


84 


DOe7 25 =S aa 
K = Kel 


B = B+C(K)*X( J) 


CONTINUE 
K = K+1 
B = BtC(K) 


WRITE(6,84) B 


FORMAT(/,2X, ‘THE SCALED MF IS :' ,2X,G15.7) 


BR = B * Z(NIND + 1) 


THE FOLLOWING ENSURES THAT THE OUTPUT STAYS IN WITHIN LIMITS. 
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XHI = 240.0 
XLO = 70.0 
BR = AMAX1(XLO, BR) 


BR 


AMIN1(XHI, BR) 


WRITE(6,85) BR 


FORMAT(/,2X, MP IS =: »2X615979) 


RETURN 


90 


C 


END 


Creare rede se devedevexe Hove ve se Fes Wee Nes FETs Fe TE Te Fe Me Ve Te TEI Ve TEM IE Ie ETE Te He Ve Fe HEE Te Ie Fe Te Fe He Fe Fee He Fe He se Fe Fe 


SUBROUTINE SUBMA(X1,X2,BR) 


Cravaeatieseaesevedevede se dede sete dete Hee ve Fede ve ae Neves de Fe se He eve He Fe Fe Fe He He VeFe Ne Hee Tee Ne He He Fe Ke Fe Fe te HN He 


C 


C 


OQ OQ © o) 


THIS SUBROUTINE PRODUCES OUTPUT 'MA' FOR THE GIVEN INPUTS. 


PAMENSTONEX(S ),CC21),2(5),XRC5) 


XR( 1) 


XR(2) 


COEFFICIENTS 


C( 
C( 
CC 
C( 
C( 
CC 


1)= 
2)= 
3)= 


4)= 


5) 


6)= 


X1 


X2 


OF THE QUADRATIC CURVE FIT. 


57 0198 

SOmy 27-0 Worl 
G6. 2529498 
QO. 1880112 
-0. 6588774 


0. 3668176 


SCALING FACTORS. 


Z(1) 
Z(2) 


es) 


36000. 0 
43.0 


13000. 0 


9] 


686 


NIND = 2 


DO 686 I = 1,NIND 
X(I) = XRCI)/2(1) 


CONTINUE 


C CONSTRUCT THE COMPLETE QUADRATIC EQUATION. 


aa 


70 


Te 


B = 0 
K = 0 
DO 70 J = 1,NIND 


DO 71 I = J,NIND 

K = K+ 
B = B+C(K)*X(J)*X(I) 
CONTINUE 


CONTINUE 


DO 72 J = 1,NIND 
K = K+l 


B = B+C(K)*X(J) 


CONTINUE 
K = Ktl 
B = B+C(K) 


WRITE(6,84) B 


FORMAT( /,2X,'THE SCALED MA IS : ,2X)Giimey 


C 
BR = B * Z(NIND + 1) 
c 
c 
C THE FOLLOWING ENSURES THAT THE OUTPUT STAYS IN WITHIN LIMITS. 
C 
C XHI = 13500.0 
é XLO = 5500.0 
c XHI = 15000.0 
a XLO = 4000.0 
c 
C BR = AMAX1(XLO, BR) 
C BR = AMIN1(XHI,BR) 
C 
C WRITE(6,85) BR 
fess) ~FORMAT(/,2X,'MA IS :',2X,G15.7) 
c 
C 
RETURN 
END 
E 


Caw eH Fever HIN Te ere Te HN Fee He HEN Fe MeN NE Te TEN TEFEN VE CN KIN HNN HEN URE KR 
SUBROUTINE SUBT2(X1,X2,BR) 

Crewe rede sede Fewer Weve Ne Weve Fede Here FeV Fe Fe WEN VEN FoI Ie He IE VE NE HE He Fe I HIER ENE TE FETTER I NI HE 

C 

C THIS SUBROUTINE PRODUCES OUTPUT 'T2' FOR THE GIVEN INPUTS. 


C 
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DIMENSION X(5),C(21),2(5),XR(5) 


XR(1) 


XR( 2) 


COEFFICIENTS 


CC 1) 


Cee) 


CU-3) 


CC 4) 


Cel) 


Cie=o) 


X1 


OF THE QUADRATIC CURVE FIT. 


SORT 1337 
2.203622 

-1. 040498 
0. 1354878 
-0. 4898891 


0. 7473461 


SCALING FACTORS. 


Zi = 36000. 0 
Ze) c= 43.0 
23) = 800. 0 
NINDo= 2 
pn DO 50051 = Ie NIND 
X(I) = XR(1I)/ZC(1) 
500 CONTINUE 
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C 


C CONSTRUCT THE COMPLETE QUADRATIC EQUATION. 


c 
Eee) 
K=0 
DO 70 J = 1,NIND 
DO 71 I = J,NIND 
K = K+l 
B = B+C(K)*X(J)*X(I) 
7 CONTINUE 
70 CONTINUE 
C 
DO 72 J = 1,NIND 
K = K+1 
B = Bt+C(K)*X(J) 
2 CONTINUE 
G 
K = K+1 
B = B+C(K) 
c 
G WRITE(6,84) B 
C 84 HORMAMG/ 2X ,enHE SCALED T2 1S ;° ,2X,G15.7) 
C 
C 
BR = B * Z(NIND + 1) 
G 
c 


C THE FOLLOWING ENSURES THAT THE OUTPUT STAYS IN WITHIN LIMITS. 
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C 

C XHI = 850.0 

C XLO = 500.0 

C XHI = 1000.0 

C XLO = 400.0 

C 

C BR = AMAX1(XLO, BR) 


C BR = AMINI(CXHI, BR) 


c WRITE(6,85) BR 


C 85 FORMATIC/ 52%, lig dione ee Gine 7) 


* 
END 
C 
Civvterededevevede dest tee vede He deve Ve Fe Fe Fe He Fe Tove HE Te Fe Vege He Me Fe Fe VEN He He He Ne He He Tee Te Fe Ne He Tesoro de Veve sede ye veN 


SUBROUTINE SUBQC(X1,X2,BR) 
Credededeveveve sede deve ve ds Fee ve te Ne ve Me eI VE TENE TE EE Fe Fe TE He IE Te FFE IE Me TEN IE VET FETE TE TE Fe FE IE I FEN Ne Te FE HE FEN 
C 


C THIS SUBROUTINE PRODUCES OUTPUT ‘QC’ FOR THE GIVEN INPUTS. 


é 
DIMENSION X(5),C(21),2(5),XR(5) 
€ 
RG) = ot 
XR(2) = X2 
C 
c 
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C 


C COEFFICIENTS OF THE QUADRATIC CURVE FIT. 


eC 
CC 
CC 
CC 
CC 
C( 
CC 
C 
C SCALING F 
C 
Aa) 
Ae> 
Z(3) 
6 
NIND 
c 
711 DO 5 


CL ) 


i=) = 92796132 
a= 202038512 
b= eo tO. 70980 
4)= 0.1464243 
5)= a0 5 7019 


6)= -0. 3884839 


ACTORS. 


36000. 0 


43.0 


F020 


ll 
NS) 


00 I = 1,NIND 


ee ely CCl) 


500 CONTINUE 


C CONSTRUCT THE COMPLETE QUADRATIC EQUATION. 


tw 
il 


Ke= 


DOr 


O J = 1,NIND 


Oy 


Gar “C2a> «C) 


O 


DO 71 I = J,NIND 
K = K+ 
B = B+C(K)*X(J)*X(1) 
71 CONTINUE 


70 CONTINUE 
DO 72-5 = 1 NaND 
K = K+l 
B= BICC Ra) 


72 CONTINUE 


Ky = Kal 


B= Bre) 


WRITE(6,84) B 


84 FORMAIC/ , 2X, THE SCALED OC 1S =", 2% isa) 


BR = B® ZCNEND + 1) 


THE FOLLOWING ENSURES THAT THE OUTPUT STAYS IN WITHIN LIMITS. 


AT = 1305.9 


XLO = 40.0 
XHI = 300.0 
XLO = 0.0 
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C BR = AMAX1(XLO, BR) 
é BR = AMIN1(XHI,BR) 
C 

C WRITE(6,85) BR 


© 85 MORMAT( /,2%, OC IS :' ,2X,G15. 7) 


C 

C 
RETURN 
END 

C 

C 

C 


Ctetede deve deve de He weve ae ve ve We se Fe ve We Ie Fe Fe ETE He He Tee We se Ve Me Ne TE Te Vee We IE ETE FETE ETE TE TEE Te IE Te Me He Te VEE Ie He 
SUBROUTINE SUBP2(X1 X42 ,X3,%4,%5 , BR) 

Crete deve dese ese Tove se ds seve Wee ve He Fe Fes Fe He He Ve Is Te Wee He Fe Ve VEN A He Fe Te TE Te We NEE TEN I Fe 76 He Te Ne Me ese Is He 

C 


C THIS SUBROUTINE PRODUCES OUTPUT 'P2' FOR THE GIVEN INPUTS. 


C 
DIMENSION X(5),C(21),Z(6),XR(5) 
C 
i 
eGo = se 
RD) 
Pe ye= x5 
MRCS) = x4 
MRS )= XS 
C 
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C COEFFICIENTS OF THE QUADRATIC CURVE FIT. 


Z(1) = 


C 
CGS) = arleaen 
Ga )= 2.899741 
Ges] “1. 2230n4 
C( 4)= -1.854803 
C( 5)= -10. 26167 
C( 6)= -0. 2169524 
C( 7)= -1.156939 
C( 8)= 2.860795 
C( 9)= 5.767990 
C(10)= -0.101891 
Cit j= none os 
C(12)= -0. 441640 
C(13)= 0. 7359644 
C(14)= -9.559825 
C(15)= 22.88968 
C(16)= 4.7794 
C(17)= -2.558953 
GC18)=) 80272274 
C(19)= 6.295503 
C(20)= -2emiav775 
C(21)= 9.380198 

C 

C 

C SCALING FACTORS, 

C 


36000. 0 


100 


C2 "A ea> 2 eae 


500 


Z(2) = 13000. 0 
Z(3) = 800. 0 
2(4) = 240. 0 
Z(5) = 2020 
Z(6) = 43.0 
NIND = 5 


DO 500 I = 1,NIND 
X(I) = XRCI)/2Z(1) 


CONTINUE 


CONSTRUCT THE COMPLETE QUADRATIC EQUATION. 


oy 


70 


DO 70 J = 1,NIND 
HOS It = J,NUND 
K = Ktl 
B = B+C(K)*X( J)*X(C1) 
CONTINUE 


CONTINUE 


DO 72 J = 1,NIND 


K = Ktl 


101 


(ys 


84 


B= BtCCK) xp 


CONTINUE 
K = K+l 
B = B+C(K) 


WRITE( 6,84) B 


FORMAT(/ ,2X, @HESSCALED P2 1S: J2x.Gioee» 


BR = B * ZCNIND + 1) 


THE FOLLOWING ENSURES THAT THE OUTPUT STAYS IN WITHIN LIMITS. 
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HT = 94S"0 

xO = 20.0 

XHI = 100.0 

Oe s0e0 

BR = AMAX1(XLO,BR) 
BR = AMIN1(XHI, BR) 


WRITE(6,85) BR 


FORMAT(/,2X,'P2 IS :',2X,G15. 7) 


RETURN 
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END 
C 
C 


Ce Mere e rede sede Wee VeVe FEE FETE He Fe FE TENE FEE FETE VE TE TE Be VE Te VE VE TE FETE TE TE TET TE HE Fe VET FE HEE HE FE VE TE HE VET HEHE HFC 


SUBROUTINE SUBT4(X1,X2,X3,X4,X5, BR) 


Creede devedededede cece Fede Te Ted TTT TVET TET DE LORE TERETE TELE TEETER LEER EE REE RERRE ERE 


C THIS SUBROUTINE PRODUCES OUTPUT 'T4' FOR THE GIVEN INPUTS. 


6 
C 
DIMENSION X(5),C(21),2(6) ,XR(5) 
C 
e 
aml) = Xl 
XR(2) = X2 
XR(3) = X3 
Res) = x4 
MRS) = X5 
C 
C COEFFICIENTS OF THE QUADRATIC CURVE FIT. 
C 
C 


C( 1)= -22.20944 
C( 2)= 10. 79398 
CS mete o9s0 
C( 4)= 86. 64350 
C( 5)= -208. 0447 


CC 6)= 1.232848 
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Cl 7 )= j-12546699 
CC 8)= -64.69914 
CC 9)= 180.0014 
C(10)=) 303647 9730 
CC ll Sod res. 
CC 12)= 420721592 
CCI3)= 16. 70037 
CQ14)= -121.1824 
CCQ15)= 183.1548 
GC 16)=  —1387 5067 
CCl] = <2117. 2744 
CC1B)=->-18202553 
CCI9}= “2275909 
CU20 y= 4229-6525 


C(21)= 73.97864 


SCALING FACTORS. 


Z(1) = 36000. 0 
Z(2) = 13000. 0 
Z(3) = 800. 0 
Z2(4) = 240. 0 
Z(5) = 20210 
Z(6) = 1800. 0 
NIND = 5 
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O 


pal 


500 


DO 500 I = 1,NIND 
X(I) = XR(1I)/2(T) 
CONTINUE 


CONSTRUCT THE COMPLETE QUADRATIC EQUATION. 


A 


70 


a2 


B= 0 
K = 0 
DO 70 J = 1,NIND 


DO 71 I = J,NIND 

K = K+tl 
B = BtC(K)*X(J)*X(1) 
CONTINUE 


CONTINUE 


DO 72 J = 1,NIND 
K = K+tl 


B = B+C(K)*X(J) 


CONTINUE 
K = Ktl 
B = B+C(K) 


WRITE( 6,84) B 
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C 84 FORMAT(/,2X,'THE SCALED T4 IS :' ,2X,G15. 7) 


BR = B * Z(NIND + 1) 


C THE FOLLOWING ENSURES THAT THE OUTPUT STAYS IN WITHIN LIMITS. 


C 
E XHI = 2000.0 

C XLO = 130000 

C XHI = 5000.0 

C Kou On0 

C 

E BR = AMAX1(XLO,BR) 
C BR = AMIN1( XHI,BR) 
c 

C WRITE(6,85) BR 


C 85 FORMAT( /#a@m@ T4 IS > ,2X,G15.7) 


C 

C 
RETURN 
END 

C 

C 


Crvvededeacae sede sede deve ae lk dese ae sea ve ae se vege Fe ese ve Hee ae Te ae Mee Te TE Fe Te TE IE TEM Fe WET Ve TET He IE IE TE IE He TE eI He 


SUBROUTINE SUBQHT(X1,X2,X3,X4,X5, BR) 


Seseskvokvoeceseseskscskstskavscsvavsicckscleskcescskslestakskatscscsksvskeskskskseskveskseskseseokdeskdeskdede ses ksksest 
C 


C 
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C THIS SUBROUTINE PRODUCES OUTPUT 'QHPT' FOR THE GIVEN INPUTS. 


C 
C 
PLMENSTON X05), C(21), 206), XRU5) 
C 
C 
XR(1) = Xl 
GZ Sexe 
XR(3) = X3 
XR(4) = k4 
XR(5) = X5 
C 


C COEFFICIENTS OF THE QUADRATIC CURVE FIT. 

c 
SC 5] Shc chink: 
C( 2)= -562. 3596 
C( 3)= -23. 65895 
C( 4)= -54.97896 
C( 5)= 98.09515 
GC Qe Oy. Bak 
G7 = 939591497 
C( 8)= 119.9962 
C( 9)= -248. 3938 
C(10)= -0,1507291 
C@ij= 17995723 
C(12)= -17. 87346 
C(13)= -40. 67739 


CQ14)= 28.27711 
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COPS )= “1902205 
C(16)= -160. 9423 
C(17)= 260. 8458 
CQ18)= 21. 40023 
CQ19)= -34. 85067 
C(20)= =21920661 


CC21)= 62. 16870 


SCALING FACTORS. 


a 


500 


Z(1) = 36000. 0 
Z2(2) = 13000. 0 
Z2(3) = 800. 0 
Z(4) = 240. 0 
Zo) 2050 
2(6) = 13020 
NIND = 5 


DO 500 I = 1,NIND 
X(T) = XRC1I)/Z(T) 


CONTINUE 
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C CONSTRUCT THE COMPLETE QUADRATIC EQUATION. 


B = 0 


K = 0 


DO 70 J 1,NIND 
DO 71 1 = J,NIND 
K = K+l 
B = B+C(K)*X(J)*X(1) 
Hal CONTINUE 


70 CONTINUE 
DO 72 J = 1,NIND 
K = K+1 
B = B+C(K)*X(J) 


ie CONTINUE 


K = Ktl 


B = B+C(K) 


C WRITE(6,84) B 


C 84  FORMAT(/,2X,'THE SCALED QHPT IS :',2X,G15. 7) 


BR = B * Z(NIND + 1) 


C THE FOLLOWING ENSURES THAT THE OUTPUT STAYS IN WITHIN LIMITS. 
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C 

C Ki le—s 13080 

C XLO = 40.0 

C XHI = 300.0 

C XO =040 

C 

E BR = AMAX1(XLO, BR) 
C BR = AMIN1(XHI,BR) 
C 

c WRITE(6,85) BR 


C 85 FORMAT(/ 52x, QHPD [sv= 52x G57) 


C 
C 
RETURN 
END 
C 
Crvdesbdcababatvedsdedeve seven ode fea ve ak aleve ve ve ve ve ve ae sede ae ve deve ve Teva ve we ede Te eI Ie Fe Te Me Fe He HE ee Ge VE 


SUBROUTINE SUBP4(X1,X2,X3,BR) 
Creve vevevedeak teaese de dk ve ve ak ae dese ak ae ae oF ve ve oF Fes OF ae aE ae IE TE TE ETE Te IE NE Te TE TEE VEE TEE TE Te Te Te ETE IE IE FE IE TE Oe 
G 


C THIS SUBROUTINE PRODUCES OUTPUT 'P4' FOR THE GIVEN INPUTS. 


C 
C 
DIMENSION X(5)),CC21)- ZG6)e RG) 
C 
C 
XR(1) = Xl 
XR(2) = K2 
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C 
C 
C 


Ca Ole Cy 


XR(3) = X3 


CORFFICIE 


CC 
CC 
C( 
CC 
CC 
CC 
CC 
CC 
CC 


CC 1 


NTS 


1) 
2)= 
3)= 
4)= 
5)= 
6)= 
7)= 
8)= 
9)= 


0)= 


OF THE QUADRATIC CURVE FIT. 


OF 19261578 

1 SeoZz6 

0. 1008366 

6. 138049E-02 
pea 29569ER-02 
so loolLale -02 
-0. 8789043 
eral? lS 
-4,834537E-02 


1.559548 


SCALING FACTORS. 


Z( 1) 
Z(2) 
Z(3) 


Z(4) 


NIND 


13000. 0 
1800. 0 
3000. 0 


20.0 


Gl 


ye 


500 


DO 500 I = 1,NIND 
X(I) = XR(I)/2Z(T) 
CONTINUE 


CONSTRUCT THE COMPLETE QUADRATIC EQUATION. 


a 


70 


v2 


84 


B= 0 
K = 0 
DO 70 J =-1,NIND 


DO 71 I = J,NIND 

K = K+l1 
B = B+C(K)*X(J)*X( 1) 
CONTINUE 


CONTINUE 


DO 72 J = 1,NIND 
K = K+l 


B = Bt+C(K)*X(J) 


CONTINUE 
K = K+l 
B = B+C(K) 


WRITE(6,84) B 


FORMAT( /,2X,' THE SCALED P4 IS). y2xGiaae 
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BR = B * Z(NIND + 1) 


C THE FOLLOWING ENSURES THAT THE OUTPUT STAYS IN WITHIN LIMITS. 


C 

e Mae e201 0 

C SALON — eed 2 

C Sai = 95080 

C XG = 0,10 

f 

c BR = AMAX1(XLO,BR) 
C BR = AMIN1(XHI, BR) 
a 

c WRITE(6,85) BR 


CG 85 Hon wipe tS ¢  .2X,G15. 7) 


C 

C 
RETURN 
END 

C 


Cretedederevede sede sede seve seve Hse Fete Ve Fe Te Ve Fe Te Fe Fe Tee Ve Fe VOTE Te TENE Te HE TE TEN Fe TEN Fe ME TE TENET NN He He HE 


SUBROUTINE SUBQFT( X1, X2,X3,BR) 


Car reresedede de deve dese Fe sede ve deeds Fede se Fe Fe Te TET He Te Te Fe Tee Me Fe Te ETE RICH RRM TR RE 
c 


C THIS SUBROUTINE PRODUCES OUTPUT 'QFT' FOR THE GIVEN INPUTS. 
C 


ii3 


C 


C 


C 


DIMENSION XC5)7,6C 21) (2Ge Ke) 


XR( 1) 
XR( 2) 
XR(3) 


COEFFICIENTS 


CC 
C( 
C( 
CC 
CC 
CC 
C( 
CC 


C( 


1)= 
2)= 
3)= 
4)= 


>) 


6)= 
7)= 
8 )= 


9)= 


C(10)= 


X1 
X2 


X3 


OF THE QUADRATIC CURVE FIT. 


2 1a 

0. 8755642 
=“Q26020919 
37 892629 

-0. 1769417 
1.446682E-02 
sige oo2s 
-7.607660 
O-2095 155 


3.747696 


SCALING FACTORS. 


Z(1) 
Z(2) 
Z(3) 


13000. 0 
1800. 0 


3000. 0 
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711 


500 


Z(4) 480.00 


It 
WD 


NIND 


DO 500 I = 1,NIND 
X(I) = XR(1)/2(1) 
CONTINUE 


CONSTRUCT THE COMPLETE QUADRATIC EQUATION. 


7B 


70 


72 


B = 0 

K = 0 

DO 70 J = 1,NIND 

DO 71 I = J,NIND 
K = Ktl 


B = B+C(K)*X(J)*X(T) 
CONTINUE 


CONTINUE 


DO 72 J = 1,NIND 
K = Ktl 
B = B+C(K)*X(J) 


CONTINUE 


K = K+l 
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B = Bt+C(K) 


C WRITE(6,84) B 


Gunes FORMAT(/,2X,'THE SCALED QFPT IS : ' ,2X,G15. 7) 


BR = B * Z(NIND + 1) 


° 


C THE FOLLOWING ENSURES THAT THE OUTPUT STAYS IN WITHIN LIMITS. 


C 

C XHI = 480.0 

C XLO = 25.0 

C 

C BR = AMAX1(XLO, BR) 
C BR = AMIN1(XHI,BR) 
C 

C WRITE(6,85) BR 


G: :5 KORMAIC/ 254-0bPt 15 425,615. 7) 


C 

C 
RETURN 
END 

C 


Crerededesevedededesedese Jove x ve ve ve ese Ne ve He Fe se ve Fe He FoF He Ve Fe Fe Fe HET Te Tee Fe TEE Te Te FETCH Fe He HEE FE Fe FE FE 


SUBROUTINE PART(A,B,MA,P2,T2,MF,NG,NS,P4,T4,MAF,WW) 


Crcdedededededededededededevesedeveves Joredevededededevescdededededededevedesededededetesededededede ve tetede Fete seve ye 


C 
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C THIS SUBROUTINE CALCULATES THE ELEMENTS OF THE 'A' AND 'B' MATRICES 


C IN THE STATE SPACE EQUATION: 


C 
C XDOT = A*X + B*U . 
C 
C 
C COMMON QC,NG,P2,QH,MA,T2,MF,P4,QF,MAF,T4,NS,QD,WW 
PMENSTON AC 3.3), BC3,2) 
REAL NG,NS,MF,MA,MAF,JG,JD,DENOM1 , DENOM2 , DENOM3 
C 
C 
C 
Gr= 0.009525 * 25* 3.14159 / 60.0 
wp— Opevese< 20283. 14159 / 60.0 
C 


C CALL SUBROUTINES TO GET PARTIAL DERIVATIVES. 


CALL 
CALL 
CALL 
CALL 
CALL 
CALL 
CALL 
CALL 
CALL 


C 


DMA(NG , P2 ,DMADNG , DMADP2) 

DT2(NG, P2,DT2DNG , DT2DP2) 

DQC(NG , P2 ,DQCDNG , DQCDP2) 
DP2(NG,MA,T2,MF,P4,DP2DNG, DP2DMF , DP2DMA , DP2DT2 , DP2DP4) 
DT4(NG,MA,T2,MF,P4,DT4DNG , DT4DMF , DT4DMA , DT4DT2 , DT4DP4) 
DQHT(NG ,MA,T2,MF,P4, DQHDNG , DOHDMF , DQHDMA , DQHDT2 , DQHDP4 ) 
DP4(MAF ,T4,NS,DP4DNS , DP4MAF , DP4DT4) 

DQFT(MAF , T4 ,NS , DQFDNS , DQFMAF , DQFDT4) 

DQD(NS , WW, DQDDNS , DQDDWW) 


weeecONEUTE THE COEFFICIENTS OF THE STATE SPACE EQUATIONS (I.E. THE 


Me? 


C ELEMENTS OF THE 'A' AND 'B' MATRICES). 


C 
J1 = DQHDMA 
J2 = DQHDMF 
J3 = DQHDT2 
J&4 = DQHDP4 
J5 = DQHDNG 
E1-= DOCBr2 
E2 = DQCDNG 
Ci = DMADP2 
G2 = DMADNG 
D1 = DT2DP2 
D2 = DT2DNG 
Al = DP4MAF 
A2 = DP4DTS 
A3 = DP4DNS 
Hl = DP2DMA 
H2 = DP2DMF 
H3 = DP2DT2 
H4 = DP2DP4 
H5 = DP2DNG 
Gl = DT4DMA 
G2 = DT4DMF 
G3 = DT4DT2 
G4 = DT4DP4 
GS = DT4DNG 
Bl = DQFMAF 
B2 = DQFDTS 
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B3 = DQFDNS 

Z1 = B2*G3*D2 + B2*G5 
Zze—s Blt B2*G2 

Z3 = Bl + B2*G1 

Z4& = B2*G3*D1 

Z5 = B2*G4 

Z6 = Z1 + 23*C2 

Z7 = Z&4 + Z3*C1 


DENOM1 = 1.0 - H1*C1 - H3*D1 


Y1 = (H5 + H1l*C2 + H3*D2) / DENOM1 
Y2 = H2 / DENOM1 
Y3 = H4 / DENOM1 


DENOM2 = 1.0 - A2*G4 


Y4 = (A2*G5 + A1l*C2 + A2*G3*D2 + A2*G1*C2) / DENOM2 
Y5 = A3 / DENOM2 
Y6 = (Al + A2*G2) / DENOM2 


Y7 = (A1*C1 + A2*G1*C1 + A2*G3*D1) / DENOM2 
DENOM3 = 1.0 - Y7*Y3 


Y8 = (Y4 + Y7*Y1) / DENOM3 


Y9 = Y5 / DENOM3 


le =aneroet Y/*Y2) -/ DENOMS 


Z8 Poeteo) Viet Zo-YS teg7*Y3"Y8 


Z9 


Bo ioe ZLIFY Sto + BS 


Gee 27 NZ ZY 10 4+ Z7*Y3*Y10 


Z10 


NSeebe fe SiC 2 + Js? 


Figlell 


2 


Sec l + J3*Die= E1 
nels) = Y11 + Y12*Y1 


ere — 2+ YI Zs Y2 


NE 


YS = 345 Yo2e ie 
G11 =" YS te oer 
Z12 =AaSs7S 

013 =a + YUS"VY 10 


FINAL FORM OF THE ‘A’ AND 'B' MATRICES. 


! NOTE ! ELEMENTS A33 AND B31 ARE NOT COMPUTED HERE BUT WERE 


DETERMINED EXPERIMENTALLY FROM GAS TURBINE TEST DATA. 


FOR ACCELERATIONS USE: 


A33 = -0.5 


B31 = 0.5 


FOR DECELERATIONS USE: 


“026 / 


i 


A33 


All = Z11/JG 


Al2 = Z12/JG 


a) IG 
A776 oD 
A22 = 29/JD 
A23 = Z10/JD 
A31 sone 
A32 ="0M0 
A33 = -10.0 
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Qa ao 


oO 


Bll = 0.0 
B12 = 0.0 


B21 = 0.0 


BZge at. On/ I) 


B31 = 10.0 


B32 = 0.0 


IF(A12. LT. 
IF(A13. LT. 
@aale ET, 
IF(A23. LT. 
OM a eae: 
iC A2 2568. 
WRITE(9,*) 
WRITE(9,*) 
WRITE(9,**) 
WRITE(9,*) 
WRITE(9 ,**) 
WRITE(9 ,**) 
WRITE(9,*) 
WRITE(9,*) 
WRITE(9,*) 
WRITE(9,*) 
WRITE(9,*) 
WRITE(9,*) 
WRITE(9,*) 
WRITE(9 ,*) 


WRITE(9,*) 


"A" AND '"B'' MATRICES FOR NG 


' 


t 


t 


.0) 
.0) 
.0) 
.0) 
m0) 
0) 


RETURN 
RETURN 
RETURN 
RETURN 
RETURN 
RETURN 


All 
Al2 
Al13 
A21 
A22 
A23 
A31 
A32 
A33 
Bae 


Beez 


All 
Al2 
Al13 
A21 
A22 
A23 
A31 
A32 
A33 
Bll 


Bae 


12] 


AND NS 


WRITE(9,*) ' B21 = ', B21 
WRITE(9,*) ' B22 = Spee 
WRITE(9,*) ' B31 = ', B31 
WRITE(9,*) ' B32 = ', B32 


WRITE(6,*) ' Ali =', All 
WRITE(6,*) ' Ad2 = Ae 
WRITE(6,*) ' Al3 Seen 
WRITE(6,*) ' AQ = “ane 
WRITE(6,*) ' A22 =" A22 


WRITE(6,*) ' A235 <= AOS 


RETURN 
END 


C 


Cetera resets te ete de Tete Fe Fee Fe ToGo ETE TE TEE LENS ETE TENE ETC TER IEEE TTT ETRE REE TERE EERE 


SUBROUTINE DMA(X1,X2,DMADNG, DMADP2 ) 
Cie te dete sede Tee deds ie Tete TIE TETS Fea ETS TEE TOTES OTE EEE EEE ER EREEREREEEEREREERRE RE 
C 


C THIS SUBROUTINE PRODUCES THE FOLLOWING PARTIAL DERIVATIVES: 


C DMA/DNG, DMA/DP2 


DIMENSION X(5),C(21),2Z(5) Xie 


XR(1) = Xl 


[22 


C2 C2 Ca Ga 


686 


XR(2) = X2 


COEFFICIENTS OF THE QUADRATIC CURVE FIT. 


C( 1)= 1.570198 
C( 2)= -0.7270151 
C( 3)= 0.2529498 
C( 4)= 0. 1880112 
C( 5)= -0. 6588774 


CC 6)= 0.3668176 


SCALING FACTORS. 


HO 36000. 0 
Z(2) = 43.0 
ee 13000. 0 
NIND = 2 


DO 686 I = 1,NIND 
M1) = RL) / 201) 
CONTINUE 


DMADNG 


2. 0*C(1)*X(1) + C(2)*X(2) + C(4) 


DMADNG = DMADNG*Z(3)/Z(1) 
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DMADP2 


C(2)*XCI) + 2. 0*€(3) x2 eee 
DMADP2 = DMADP2*2Z(3)/Z(2) 


C 

C 

C 
RETURN 
END 

C 


Crererereve dese sevevtede dene ae ae ae vee seve se seve ve vevede ve se Weve eve Te vee Te Ve ye HE VET Me Fe VE TE HE Fe IE Ve se He Ve TEE Te Ie VE HE 
SUBROUTINE DI2(Cxi x2. DiZENG.UIZbe2 
Credsasae serene de vesede ae ede dese ve vege ve Fe ve Ve Te Fe Ne Ve TEV VEC Ve Ie VEE TE Fe Ie NE Te VEN Fe NE WEE Ve TENE Te Ve Te Ve se He FEN Ve NE 


C 


C THIS SUBROUTINE PRODUCES THE FOLLOWING PARTIAL DERIVATIVES: 


C DIZ/DNG> Di27DPZ 


DIMENS TON 405 }5CC21) 32052 aks) 


XRCID al 


X2 


XR( 2) 


C COEFFICIENTS OF THE QUADRATIC CURVE FIT. 


CQe1) =" -0n or 1597 
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CO 2)=erZ. 203628 

Cc 3)= -1.040498 
CC 4)= 0.1354878 
CG) = =054695891 


CC 6)= 0. 7473461 


C SCALING FACTORS. 


C 
Z(1) = 36000. 0 
Z(2) = 43.0 
Z(3) = 800. 0 

C 
NIND = 2 

C 

(at DO 500 I = 1,NIND 


m@bjo= ARC1)/2U1) 


500 CONTINUE 


C 

G 

C 
DT2DNG = 2.0*C(1)*X(1) + C(2)*X(2) + C{4) 
DT2DNG = DT2DNG*Z(3)/Z(1) 

C 
DT2DP2 = C(2)*X(1) + 2. 0*C(3)*X(2) + C(5) 
DT2DP2 = DT2DP2*2(3)/2(2) 

c 

C 


RETURN 


END 
C 


Creve de deve dee ve He Fede He te ese de Ted Fe He He eH Te Ie NE TIEN Ie Te ETE TE NE TE He Te ETCH Te HI Te HET TE EK KEK 


SUBROUTINE DQC(X1,X2,DQCDNG , DQCDP2) 


Crete Here de deve dete ve Te vous Fete Te Te Fe Nee Hee Ne eH Fe HITE IEA IEICE TET HT ICN HE TIEN I III RI 


C 
C 
C THIS SUBROUTINE PRODUCES THE FOLLOWING PARTIAL DERIVATIVES: 
C 
C DQC/DNG, DQC/DP2 
C 
C 
DIMENSTON X(5) ,€(C21),.2( 5), XKGa) 
G 
XR(1) = Xl 
XR(2) = X2 
C 
C 
C 


GC COEFFICIENTS OF THE QUADRATIC CURVE FIT. 


C 
CC 1) = 9-97 263 2 
Cl 2)= 520703512 
CC 3)= -10. 70980 
CC 4)= 0.1464243 
CC 5)= eGo? og 
CC 6)= -0. 3884839 

C 
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C SCALING FACTORS. 


C 


Z(1) = “6000. 0 
Z(2) = 43.0 
Z(3) = T3070 


Hi 
ho 


NIND 


poll DO 500 I = 1,NIND 


X(I) = XRC(I)/ZCI1) 


300 CONTINUE 


2.0*C(1)*X(1) + C(2)*X(2) + C4) 


DQCDNG 


DQCDNG = DQCDNG*Z(3)/2Z(1) 


DQCDP2 Gez = xC1) + 2.0*C(3)*X(2) + C5) 


DQCDP2 = DQCDP2*Z(3)/Z(2) 


RETURN 


END 


C 
C 


Cvededede dete de dete de dete fede fede de Fe Bede Fe FC PT OLE EVEL TE PERE TERE PEALE TERE EE RTE REPT ER EEE EEE ERER 


SUBROUTINE DP2(X1,X2,X3,X4,X5 ,DP2DNG ,DP2DMF ,DP2DMA , DP2DT2 , DP2DP4 ) 


Cretedededevededededstevedede ede ds dete Fe vedo tededetetetedededetecesetedesetedededscededesedece le dedede dete lekR FRIAR AAG 08 


C 


C THIS SUBROUTINE PRODUCES THE FOLLOWING PARTIAL DERIVATIVES: 


C 
C DP2/DNG, DP2/DMF 
C 
C 
DIMENSION X(5),C(21),2(6),XR(5) 
C 
C 
RRC = od 
eee? 
Mes — oe 
XR(4) = X4 
XR(5) = X5 
E 


C COEFFICIENTS OF THE QUADRATIC CURVE FIT. 
c 

C( 1)= 4.17287 

C( 2)= -2.839741 


C( 


WW 
/ 
ll 


le 2230146 
CC 4)= -1. 854803 
CC 5)= -10. 26167 
C( 26) = 5-022 1169524 
COM = = els 6937 
CC -8)= -2, 860725 
C{ 9)=5 767990 
C( 10) =se=0101891 


CCl1)= 0.2918 
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C2 Gl ea 


C(12)= -0. 441640 
C(13)= 0. 7359644 
C(14)= -9.559825 
C(15)= 22. 88968 
C(16)= 4.7794 
C(17)= -2.558953 
C(18)= 0.272224 
C(19)= 6.295503 
C(20)= -28.57775 


C(21)= 9.380198 


SCALING FACTORS. 


yi 


500 


Z(1) = 36000. 0 
Z(2) = 13000.0 
Z(3) = 800. 0 
Z(4) = 240.0 
Z(5) = 20. 0 
Z(6) = 43.0 
NIND = 5 


DO 500 I = 1,NIND 
X(T) = XR(1)/2(T) 
CONTINUE 


C 


C 


DP2DNG 


DP2DNG 


DE2CoHE 


DP2DMF 


DP2DMA 


DP2DMA 


DE 2ZDT2 


DP 2012 


DP2DP4 


DP2DP4 


RETURN 


END 


2*CC1)*X(1) + CC2)*X(2) + CCS) 3X03) Gee ae 
+ COS )EKC e+ CC cy) 


DP2DNG*Z(6)/Z(1) 


C(4)*X(1) + CC8)*X(2) + C(11)*X(3) + 2*C(13)*X(4) 
+ C(14)*X(5) + C(19) 


DP2DMF*Z(6)/Z(4) 


C(2)*X(1) + C(7)*X(3) + C(8)*X(4) + 2*C(6)*X(2) 
#7C09)7x05), 1 CO) 


DP2DMA*Z2(6)/Z(2) 


C(3)*X(1) + CC7)*X(2) + CC11)*X(4) + 2*C(10)*X(3) 
+ C(12)*xX(5) + C(18) 


DP 2012206) 7/203) 


C(C5)*XC1) + CC9)*X(2) + CC12)*X( Ge 2* Caras em 
+ €(14)*X(4) + C(20) 


DF 2DP4*2(6)/Z(5) 


Cretededevedededededededede ce eve dete Feds Dede Bee TE TG ETE FCTE TT BETTE TEER TERETE BER RETR ETE TREE EE EERRE CE 
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SUBROUTINE DT4(X1,X2,X3,X4,X5,DT4DNG, DT4DMF , DT4DMA , DT4DT2 , DT4DP4 ) 


Crededesesesededese secede decease Fev eae Fe ToT NE TENE Vee Te Ie IE Te VE TE TE Te IE TE TENE IE Te HEHE Te TC LETC NE Te ESTE TE Fe TeV TENE TC TERE TEE Te Ie Ge Tee Se 


E 
C 
C THIS SUBROUTINE PRODUCES THE FOLLOWING PARTIAL DERIVATIVES: 
C 
C DT4/DNG, DT4/DMF 
C 
C 
C 
DIMENSION X(5),C(21),2(6),XR(5) 
C 
C 
eG el 
MRC) a= X2 
XR(3) = X3 
XR(4) = X4 
XR(5) = X5 
E 


C COEFFICIENTS OF THE QUADRATIC CURVE FIT. 


CC 1)= -22. 20944 
CC 2)= 10. 79398 
C( 3)= 21.99301 
CC 4)= 86.64350 
CC 5)= -208. 0447 


CC 1. 232848 


ON 
ww 
il 
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CC 7)= 
CC 8)= 
CC 9)= 
C(10)= 
C(11)= 
C(12)= 
C(13)= 
C(14)= 
C(15)= 
C(16)= 
C(17)= 
C(18)= 
C(19)= 
C(20)= 


C(21)= 


-12. 46899 
-64. 69914 
180. 0014 
-0. 6479730 
=1iOte Is 
20m Zio 
16. 70037 
-121. 1824 
183. 1548 
130.9007 
= ie a 
O20 95 95 
12.759 89 
=229 4355 


73. 97864 


SCALING FACTORS. 


Z(1) 
TD) 
73) 
Z(4) 
Z(5) 
Z(6) 


NIND 


ll 


36000. 0 
13000. 0 
800. 0 
240.0 
20.0 


1800. 0 


[32 


vA 1 


500 


DO 500 I = 1,NIND 
CD Se MR AC 


CONTINUE 


DT4DNG = 2%*C(1)*X(1) + C(2)*X(2) + C(3)*X(3) + C(4)*X(4) 
+ C(5)*X(5) + C(16) 


DT4DNG = DT4DNG*Z(6)/2Z(1) 


DT4DMF = C(4)*X(1) + C(8)*X(2) + C(11)*X(3) + 2*C(13)*X(4) 
+ 0(14)*X(5) + C(19) 


DT4DMF = DT4DMF*2Z(6)/Z(4) 


DT4DMA = C(2)*X(1) + C(7)*X(3) + C(8)*X(4) + 2*C(6)*X(2) 
Smee G5) +) CC 17) 


DT4DMA = DT4DMA*Z(6)/Z(2) 


DT4DT2 = C(3)*X(1) + CC7)*X(2) + C(11)*X(4) + 2%C(10)*X(3) 
+ €(12)*X(5) + C(18) 


DT4DT2 = DT4DT2*Z(6)/Z(3) 


Dieweam—sGls)*xX( 1) + CC9)*XC2) + CC12)*XC3) + 2*C(15)*X(5) 


+ €(14)*X(4) + C(20) 


DT4DP4 = DT4DP4*Z(6)/Z(5) 


53 


RETURN 
END 


PERETTI TET TERETE TCT TCP TER TE BETTER ERE ECE RRR REE RERIERERERER ER ERERERERS SM 15510 


SUBROUTINE 
X2,X3,X4,X5 , DQHDNG , DQHDMF , DQHDMA , DQHDT2 , DQHDP4)SSM15520 


ede desea de ae ae de de ae ae de ae aie aie aie ae aie ak ae ale se ae aie ae ak de ae ae ae oe ae aie aie ae aie ale ale aie aie ake ae ae He Te Fe He ae ae ae ak ae He ak He ae He ae Ie eee SSM 15530 
c 


C 


C THIS SUBROUTINE PRODUCES THE FOLLOWING PARTIAL DERIVATIVES: 


C 
c DQH/DNG, DQH/DMF, DQH/DMA, DQH/DT2, DQH/DP4 
G 
c 
e 
DIMENSION X(5),C(21),2(6),XR(5) 
c 
E 
KCl 
MRG2 =e 
XR(3) = X3 
MCG = x4 
KRG 00 
C 


C COEFFICIENTS OF THE QUADRATIC CURVE FIT. 


CC 1)= 343.8178 
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CC 
CC 
CC 
C( 
CC 
CC 
CC 
CC 
C(1 
C(1 
C(1 
C(1 


CCl 


0)= 
1)= 
2)= 
3)= 
4)= 


C(15)= 


CCl 


C(1 


6)= 


7)= 


C(18)= 


C(19)= 


C(20)= 


¢(21)= 


Ol. 3590 
= Sin lepeneye ls. 
-54. 97896 
GSO S SS 
217. 8508 
3, 591497 
gE 996? 
-248. 3938 
Oyo O7 291 
WeI>/22 
-17. 87346 
-40. 67739 
Te PT TAU 
19 Oe 2205 
-160. 9423 
260. 8458 
Zea o0 23 
-34. 85067 
ao Ooo t 


62. 16870 


SCALING FACTORS. 


Gi) 
Z@2) 
Z(3) 
Z(4) 


36000. 0 
13000. 0 
800. 0 


240.0 
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Z2(5) = 20.0 
Z2(6) = 15020 
NIND = 5 


DO 500 1 = 1 ININD 
ACE) SAR aa) 


CONTINUE 


DQHDNG = 2*C(1)*X(1) + C(2)*X(2) + C(3)*ACS) + CCA es 
tuGG5)% Go) Colo) 


DQHDNG = DQHDNG*Z(6)/Z(1) 


DOHDMF = C(4)*X(1) + C(8)*X(2) + C(11)*X(3) + 2*C(13)*X(4) 
+ C(14)*X(5) + C019} 


DQHDMF = DQHDMF*Z(6)/Z(4) 


DQHDMA = C(2)*X(1) + C(7)*X(3) + C(8)*X(4) + 2*C(6)*X(2) 
+ C(9)*X(5) + C(17) 


DQHDMA = DQHDMA*Z(6)/Z(2) 


DQHDT2 = C(3)*X(1) + C(7)*X(2) + €(11)*X(4) + 2*CCIG)a es 


+ €( 12)*X(5) 4 Cis) 


DQHDT2 = DQHDT2%*2(6)/2(3) 
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DQHDP4 ey ye J ee 12 )4K03) + 2*C(15)*xX(5) 


1 PCC x4) Pac C20) 


DQHDP4 = DQHDP4*Z(6)/Z(5) 


RETURN 
END 
C 


Crevededeve sede fede deve fede deve FER TCT TS ETE TET TE TCR TOTES OTE TE TERETE TET TERS DTS BRE TER TTR TTE 


SUBROUTINE DP4(X1,X2,X3,DP4DNS , DP4MAF , DP4DT4 ) 


THIS SUBROUTINE PFODUCES THE FOLLOWING PARTIAL DERIVATIVES: 


DP4/DNS 


oO o 2 C2 2 © © © 


DIMENSION X(5),C(€21),2(6) ,XR(C5) 


XR(1) 


X1 


XR( 2) X2 


XR(3) = X3 


C COEFFICIENTS OF THE QUADRATIC CURVE FIT. 


7 


Con) 
c( 2) 
c( 3) 
c( 4) 
c( 5) 
C( 6) 
C( 7) 
Bere) 


CC 9) 


C(10)= 


ce 


1926175 

yl os2e 

. 1008366 

. 138049E-02 
429369 =02 
5. 13614 lhe o2 
0. 8789043 


Le Sa 


-4.6345375-02 


1 


- 559548 


SCALING FACTORS. 


500 


Z(1) 


Z(2) 


Z2(3) 
Z(4) 


NIND 


DO 500 


3 


= 


13000. 0 
1800. 0 
000. 0 


2028) 


1,NIND 


ACT) 9 = Re 2a 


CONTINUE 
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DP4DNS = C(3)*XC1) + CC5)*XC2) + 2*C(6)*X(3) + CC9) 
DP4DNS = DP4DNS*2Z(4)/Z(3) 


DP4MAF = C(2)*X(2) + C(3)*X(3) + 2*C(1)*X(1) + C(7) 


DP4MAF = DP4MAF*Z(4)/Z(1) 


DP4DT4 = C(2)*X(1) + CC5)*XC3) + 2*C(4)*XC2) + C(8) 


DP4DT4 = DP4DT4*Z(4)/Z2(2) 


RETURN 
END 

C 

Crete rete secede Feet eT PETS PETE TELE TE TE TE TEE TE ETE TS TET TSEC TS TSE TC TET TE TS TTT TERETE TT ATER BREE RE 
SUBROUTINE DQFT(X1,X2,X3,DQFDNS , DQFMAF , DOFDT4) 

Cede cede te vedede feted tetera sede Be Ne Bose TE TE Te OPES OEE ETE TERE TOPE PELE TELE SESE LEEPER ELE LEER RE 


C 


C THIS SUBROUTINE PRODUCES THE FOLLOWING PARTIAL DERIVATIVES: 


C DOF/DNS, DQF/DMAF, DQF/DT4 


©o2 2 © 


DIMENSION X(5),C(21),2(6) ,XR(5) 
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XR( 1) 


RREZD 


XR(3) = 


COEFFICIENTS 


C( 1)= 
C( 2)= 
C( 3)= 
Oh 
c( 5)= 
C( 6)= 
C( 7)= 
C( 8)= 
C( 9)= 


C(10)= 


X1 
X2 


X3 


OF THE QUADRATIC CURVE FIT. 


Ze 1 D248 77 

0. 8755642 
“026026919 
O76 92829 

=0, 769017 
1. 446682E-02 
=e O32. 

=7 007 660 
022095135 


3. 747696 


SCALING FACTORS. 


Z(1) 


7 


Z(3) 


Z(4) 


1300070 
1800. 0 
3000. 0 


480. 00 
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NIND = 3 


yAl DO 500 I = 1,NIND 
X(I) = XRCI)/ZCI) 


500 CONTINUE 


C 

C 

E 
DQFDNS = C(3)*X(1) + C(5)*X(2) + 2*C(6)*X(3) + C(9) 
DQFDNS = DQFDNS*Z(4)/Z(3) 

C 
DQFMAF = C(2)*X(2) + C(3)*X(3) + 2*C(1)*X(1) + C(7) 
DQFMAF = DQFMAF*Z(4)/Z(1) 

c 
DQFDT4 = C(2)*X(1) + C(5)*X(3) + 2*C(4)*X(2) + C(8) 
DQFDT4 = DQFDT4*Z(4)/2Z(2) 

c 

c 

C 
RETURN 
END 

c 


Cite te devedese te ete Rede TTT RT TT TTT TET TTT TEE RTE ETRE REREERETESERTTEERR ICT 
SUBROUTINE DQD(X1,X2,DQDDNS , DQDDWW ) 

Crete deve decode dese sede ede Fe eR TD TEL TT EERE RETR ETRE ER RRR EREREREREERRRERE 

C 


C 
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C THIS SUBROUTINE PRODUCES THE FOLLOWING PARTIAL DERIVATIVES: 
C 


C DQD/DNS, DQD/DWW 


C SCALING FACTORS 


C 
XQD = 480. 
XNS = 3000. 
XWW = 49. 
C 
Clo =-=20e80 
C2 = 4. OE-6 
Go = 121929685 
C 
BDQBDNS = 2*xX1C2 9203 (x25 3) 
C 
DQDDWW = 1. 3*C3*X1*X1*( X2**0, 3) 
C 
C DQDDNS = DQDDNS*XNS/XQD 
C DQDDWW = DQDDWW*XwWW/XQD 
C 
C 
C 
RETURN 
END 
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APPENDIX E 


TITLE GT DYNAMIC PROGRAM 
v 


we 
ae 


Jevededededededededesevesesisededevedevedekstededededededevesesedeektededesetescddededesesecesesededidetededededededevete Rae ete RR 


ve =e 
aC ve 
* BOEING MODEL 502-6A GAS TURBINE * 
a ve 
* DYNAMIC COMPUTER SIMULATION * 
ae w¥ 
ve MODIFIED BY * 
# V.A. STAMMETTI * 
FROM A ROUTINE BY 

V.J. HERDA a 
7 ec 
* THIS PROGRAM SIMULATES THE DYNAMIC RESPONSE OF THE NPS * 
* BOEING GAS TURBINE TEST FACILITY USING A MULTILPLE * 
* LINEARIZATION TECHNIQUE. * 
v ve 


Wistevedededetedetetedesedevesedededeke Fever teseve tee seve sede sete te fede Fekete FoR RNR RRR RRR RR TERRE 


C 
C 


PARAM JG=0.009525, JD=0.6738, PI=3.14159, T = 2.0 , TW=2.0 , TW1=1.0 
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xk 


* 


a 


* 


rig 


C 


C 


THE FOLLOWING VALUES LISTED ON THE FUNCTION CARD ARE FOR FUEL FLOW, 
GAS GENERATOR SPEED, AND DYNO SPEED AS A FUNCTION OF TIME. 
THESE VALUES WERE OBTAINED FROM STRIP CHART RECORDS AND ARE ENTERED 
IN THE FORM (E.G. FUEL FLOW) ...TIME(SEC), EUEESEEO} eee 


THIS SET 1S FOR EXPERIMENTA RUNG? aie 


THIS SEP IS FOR BAPERTHENTAL RUNS a. 


AFGEN NGDATA= 0.0,34900.0, 5.0,34900.0, 10.0,34900.0, 15.0,34900. OF se 


20.0,34900.0, 25.0,34900.0, 30.0,34900.20,, 35.0 34900a0 


AFGEN NSDATA= 0.0,570.0, 1.0,570.0, 3.0,57070) 5. Oy) 07 eee 


6.0,653.0, 7.0,708.4, 8.0,791.4, 9.0,846.7, 10.0,929.8, ... 
11.0,1040.5, 12.0,1178.8, 13:@903449, 14.0,1566,0 98 
15.0;,1787.7, 16.0,2009.1, 174@2230,5 Maile. 0, 245 aeeR 
19.0,2590. 2, 20.0,2673. 3, 2is0n27560. 322, 0, 2899 ee 
23.0,2894.65, 24.0,2950.0, 25.0,2950.0, 26.0) 2950 "0mm 


35.0,2950. 0 


AFGEN MFDATA= 0.0,193.5, 4.0,193.5, 5.0,195.6, 6.0,197.7, ... 


7.0,197.7, 80,197. 7, 9:0,199: 75, 110) Ge 100uSIee 
11.0,199.75, 12.0,201.8, 13.0,201,8, 1aeomegs soe 
15.0,203.9, 16.0,206.0, 20/0.206.0,) 2500 —20c"0m——nE 


30.0;20680). 35. ON2Z0630) 


AFGEN QDDATA= 0.0,450.0, 4.0,450.0, 5.0,445.1, 6.0,445.1, .. 


7.0,435.4, 8.0,425.6, 9.0,415.8, 10,0,400,1 0 


11.0,391.1, 12.0,376.8, 13.0,36@2 1, 1400nos2 0 eee 
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mmmecigs S elGuoe2es. 9, 17.0.274.3, 18.0,264.5, ... 


oreo or C00. 254. 68 21, 0,245.0, 25.0,245.0, 


30.0,245.0, 35.0,245.0 


¥ 
¥ 


INITIAL 


a 


* ESTABLISH INITIAL CONDITIONS. 


NGI=34900. 0 
NSI=570.0 
MFI = 193.5 


QDI = 450.0 


ve 


al. 


* ESTABLISH FINAL CONDITIONS 
ve 

NGF=34900. 0 
NSF=2950.0 

QDF = 206. 00 


MFF = 245.0 


* SET INITIAL STATE PERTURBATION TO ZERO 


DNG = 0.0 

DNS = 0.0 

DEF =O. 0 
¥ 

EQ = MFI 
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* 


NGB = 26000. 0 
NSB = 1750.0 
* 
* 
DYNAMIC 
¥ 
* DATA CURVE FORMULATION 
* 
NGD = AFGEN(NGDATA ,TIME) 
NSD = AFGEN(NSDATA , TIME) 
MFD = AFGEN(MFDATA, TIME) 
QDD = AFGEN(QDDATA, TINE) 
* STATE SPACE LINEAR MODEL FORMULATION 


All = -1. O*EXP(-6. 9929E-01*(NSL/NSB) + 5. 585157007 (NGE Veber 


= 3.24 53b 500) 


Al2 = -1. OWEXP(-2. 8415E+00*(NSL/NSB) + 7.9978E+00*(NGE/ NGE ee 


-4. 4662E+00) 


* A13 = (2. 1503E+03)*((NSL/NSB)**2.0) ... 

ve + (-5.9752E+0:)*(NSL/NSB)*(NGL/NGB) . 

* + (-5, 0697E+03)*((NGL/NGB)**2.0) + (1.3101E+03)%*(NSL/NSB) ... 
ve + (1. 8551E+04)*(NGL/NGB) - 9.1460E+03 


A13 = EXP(-0.45788%*(NSL/NSB) + 1.189%*(NGL/NGB) +6. 8305) 


ve 
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a 


ds 


x 


ay 


Ye 


ie 


vc 


Pee eb o tone ies CNSL/NSE 2.0) ... 


+ (-9.5351E-01)%*(NSL/NSB)*(NGL/NGB) ... 


+ (1.5745E+00)*( C(NGL/NGB)**2.0) + (5. 1810E-01)*(NSL/NSB) ... 


Pa lo2 52h 00) sCNGE/NGE) + 4.6015E-01 


A22 = (5. 6875E-02)*(NSL/NSB) + (-1.3166E+00)*(NGL/NGB) ... 
+o. JBO2E 0m 
A23 = (-1. 7434E+01)*((NSL/NSB)**2.0) ... 


+ (3.5345E+01)*(NSL/NSB)*(NGL/NGB) ... 


+ (4.5787+01)*((NGL/NGB)**2.0) + (1.0894E+01)*(NSL/NSB) ... 


+ (-2,0510E+02)*(NGL/NGB) + 1.5265E+02 


A23 = EXP(0.92011*(NSL/NSB) -4. 2549*(NGL/NGB) + 6. 2290) 
A31 = 0.0 

Bee — i050 

A33 = -10.0 

B11 = 0.0 

pie = 0.0 

B21 = 0.0 

E22 = -14. 1723156 
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B32 = 0.0 
* 
A= 1.0 
B = -2.0*(Al1 + A22) 


C Al1*A22 - A21*A12 


IMCHK = B¥*2.0 - 4. O*A*C 


IFC IMCHK. LT. 0. 0)° iMCHK=—= vee 


LAMDA1 (=B ae SORT CINCHER 7725707 A 


LAMDA2 = (-B - SQRTCIMCHK))/2.0/A 


DERIVATIVE 


NOSORT 


*« “COMPUTE INPUTS TO THE MULTIPLE LINEARIZATION ODEs. 


* RUN #1 


Dit =o Deis 
DQD = QDD - QDI 


DNGDOT Al1*DNG + A12*DNS + A13*DE 


DNSDOT = A21*DNG + A22*DNS + A23*DE + B22*DQD 


DEDOT = A33*DE + B31*DMF 


* DYNAMIC EQUATIONS FOR MULTIPLE LINEARIZATION MODEL. 
DNG=INTGRL( 0. 0 , DNGDOT) 


DNS=INTGRL(0. 0, DNSDOT) 
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DE =INTGRL(0. 0,DEDOT) 


NGL = NGI + DNG 


NSL = NSI + DNS 


SORT 


* THE STATEMENTS IN THE PREVIOUS (DERIVATIVE) SECTION YIELD VALUES 


* OF 'NG', AND 'NS' AS CALCULATED MULTIPLE LINEARIZATION MODEL. 


* MODELS. THE STATEMENTS BELOW ARE THE DSL STATEMENTS THAT SPECIFY 


% THE VARIOUS SIMULATION SETTINGS. 
ve 
TERMINAL 


METHOD RK5 


Pee oEBRR NS = 1.E-6, NG = 1.E-6, DNS 


Meeps RR NS = 1.E-5, NG = 1.E-5, DNS = 1.E-5, DNG 


CONTROL FINTIM=35.0,DELT=0. 001 
MeermiNt .5,NS,QFPT,T4 


C PRINT 0.5,NS,NSL,NSD,NG,NGL,NGD,E,WWM 


PRINT . 35,NGD,NGL,NSD,NSL 

C PRINT .5,DNG,DNS,DE 

C PRINT .5,P2GI,P4GI,P2,P4,NSD,NS,NSL,NG,NGL 
C PRINT .5,QD,QDM,QHPT,QFPT,QC,E,EF,MFM 

* SAVE 1.0,MFM,NG,NGD,NS,NSD 

SAVE .05,MFD,NGD,NSD,NGL,NSL,E,EF,QDD 


ve 


* RUN #1 


ve 


1.E-6, DNG 


EO 


15 


* GRAPH (DE=TEK618) TIME(LO=0.0,SC=0. 2,TI=.50,NI=10,UN=SEC) ... 
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af 
at 


x 


ve 


ac 


+ 


RUN #4 


GRAPH (DE=TEK618) 


RUN #4 


* GRAPH (DE=TEK618) 


RUN #1 


NS(LO=900 ,SC=100, TI=1. 6,NI=5,UN=RPM) ... 
NSD( LO=900 , SC=100 , TI=1. 6 ,NI=5,UN=RPM) ... 
NSL( LO=900 ,SC=100,TI=1. 6,NI=5 ,UN=RPM) ... 
FF(LO=100,SC=10. , TI=2. ,NI=4 , UN='LB/HR’ ) 


MFM( LO=100,SC=10. ,TI=2. ,NI=4,UN='LB/HR' ) 


NG( LO=24000 , SC=1000 , TI=1. 1428,NI=7,UN=RPM) ... 
NGD( LO=24000 , SC=1000, TI=1. 1428,NI=7,UN=RPM) .. 


NGF( LO=24000,SC=1000,TI=1. 1428,NI=7,UN=RPM) . 


TIME( LO=0. 0,SC=0. 2, TI=. 50,NI=10,UN=SEC) ... 


NG( LO=24000 ,SC=1000, TI=1. 33,NI=6,UN=RPM) ... 


NGD( LO=24000 ,SC=1000 , TI=1. 33,NI=6,UN=RPM) ... 


NS( LO=700 ,SC=100 ,TI=1. 6,NI=5,UN=RPM) ... 
NSD( LO=700 ,SC=100, TI=1. 6,NI=5,UN=RPM) ... 


MFM(LO=100,SC=10. ,TI=1. 6,Ni=5,UN= Levi ) 


TIME( LO=0. 0,SC=0. 2,TI=. 50,NI=10,UN=SEC) .. 


NG( LO=21000 ,SC=1000 , TI=1. 33,NI=6,UN=RPM) ... 


NGD( LO=21000 ,SC=1000 , TI=1. 33,NI=6,UN=RPM) ... 


NS( LO=900 ,SC=100 , TI=2.0,NI=4,UN=RPM) ... 
NSD( LO=900 ,SC=100,TI=2.0,NI=4,UN=RPM) ... 


MFM( LO=80 ,SC=10. ,TI=2. 0,NI=4 ,UN= LB/HR ) 


THIS IS FOR THESIS PRESENTATION FIGURES 
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¥ 

GRAPH (G2, DE=TEK618) TIME(LO=0. 0,NI=15,TI=.50,UN=SEC) ... 
NSD( LO=250, SC=100 ,TI=. 25,NI=36,UN=RPM) ... 

* NS(10=500,5C=100,TI=. 25 ,NI=36,UN=SEC) ... 
NSL( LO=250 , SC=100,TI=. 25,NI=36 , UN=RPM) 

* BRG@LO=50) ,SC=10sal—2aNi=4.UN= LB/HRe)) on 

* MFM( LO=80. ,SC=10. , TI=2. ,NI=4,UN='LB/HR' ) 


GRAPH (Gl, DE=TEK618) TIME(LO=0.0,TI=.50,NI=15,UN=SEC) ... 


NGD( LO=20000 , SC=1000 ,TI=. 25,NI=36,UN=RPM) ... 


* GC LO=30000 ,SC=1000 ,TI=. 25,NI=36,UN=RPM) ... 
NGL( LO=20000 ,SC=1000, TI=. 25 ,NI=36 , UN=RPM) 


GRAPH (G3, DE=TEK618) TIME(LO=0.0,TI=.50,NI=15,UN='LB/HR') ... 


¥ E(LO=180.0,SC=10 .,TI=.5 ,NI=14,UN='LB/HR’) ... 
* Et LO—lclmmsc—100il—.5 -NI—14UN='LB/HR ) ... 
¥ MEDCLO—loOnesG—10 = 1i=. 5° .NI-14,UN= LB/HR’ ) 


* GRAPH (G2, DE=TEK618) TIME, NSD, NS, NSF 
* GRAPH (G1, DE=TEK618) TIME, NGD, NG, NGF 
LABEL (G1) GAS GENERATOR SPEED 

LABEL (G2) DYNAMOMETER SPEED 

C 

END 

STOP 


C FORTRAN 
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APPENDIX F 


COMPARISON OF "A" MATRIX COEFFICIENTS 


Tables 4-9 compare the individual state space "A" matrix 
coefficients used in the dynamic computer simulation to 
those obtained from the Smoothed Dynamic Transition Map. 
The equation used by the dynamic simulation to calculate a 
particular coefficient 1s specified in each table by note 3. 
The scaling factors NSB = 1750 rpm and NGB = 26000 rpm were 
used to scale the variables NSL and NGL for a two variable 


linear regression in the forms shown. 


Te 


NG = 37000 


E> OCe 


32000 


30000 


25000 


23000 


22000 


20000 


17000 


15000 


NOTES: 


TABLE 4 


AieMArREx CORRE LEIENT FCORRELATION DATA 


RPM 


NS = 500 RPM 1000 


0): 


=OO8 


30; 
(-40. 


=20% 
C=25e 


Or 


=o. 
(a2. 


20. 


Zz) 


71 


83 
00) 


06 
00) 


86 


00) 


. 46 


60 
00) 


. 34 


a23 


80 


=oo iO) 


-48. 08 


-25,25 
(-25.00) 


-16. 43 
(-20. 00) 


-5. 62 
(-5. 00) 


p05 


a2 95 
(G=2..20) 


sage 


a0 


=0. 65 


1500 


O10: 


Sore 


=20, 


Cellier 


=13. 


(-14. 


-4, 
(-4. 


aie 
C2 


=O; 
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of 


67 
00) 


45 
DO.) 


0 
50) 


9 


41 
20) 


od 


Oo 


54 


2000 


=49), 


Sie 


mil ©: 
Co. 


= 11. 
(C= 10: 


=o); 
(-4, 


Sal. 
CZ, 


a0: 


54 


34 


93 
00) 


02 
00) 


a 
50) 


»45 


98 
10) 


oo 


. 68 


a4 


2500 


-40. 56 


-26. 40 


-13. 86 
(-12. 00) 


=o. 02 
i= 9/700 ) 


moO 8 
(-4.50) 


= 2201 


-1. 62 
(-2. 00) 


one 


=O 255 


=O. 


VALUES WITHOUT PARENTHESES ARE SIMULATION VALUES. 
VALUES IN PARENTHESES ARE VALUES FROM THE SMOOTHED 


DYNAMIC TRANSITION MAP. 


THE TWO VARIALBLE LINEARLY REGRESSED EQUATION FOR 
MATRIX COEFFICIENT All IS: 


All = -1. O*EXP(-6.9929E-01*(NSL/NSB) + 
5.5831*(NGL/NGB) - 3.2433) 


153 


TABLE $5 


Al2 MATRIX COEFFICIENT CORRELATION DATA 


NS = 500 RPM 1000 1500 2000 2500 

NG = 37000 RPM -447. 39 =19Gm05 =68721 =SOnele -l7 ae 
35000 “2G 82 =v 56 -47. 68 a2 ale -9. 40 
32000 =910. 09 -42.67 =18 39> -8.41 -3 974 


(=70' 00) -357. 00) (=5.450) (-1280o) ( -295) 


30000 So 94 =23..06 =10.°24 “4.55 -2702 
(=50. 00) Mea -25 70 at 1220) (-4. 60) ( -27003 


25000 =H 16 -4.95 Le, =O5 26 -0. 43 
(275,100) C-7 3003) Gores 3) C= 0). 90) (-0°5c5 


23000 =6. 03 ae OG ge ke, “U7 5 -0: 25 


22000 -4.43 le OY, =). 37 =O7 34 -O0 
C=5.00) L=0- 60) C050) (-0. 40) (-0, 205 


20000 =2.:39 = 206 =0..47 =02 21 -O709 
17000 a0R.2>5 0.42 0 0,706 -0.04 
15000 20-51 205.25 =“O210 =(. 05 -0.-02 


NOTES: 1). VALUES WITHOUT PARENTHESES ARE SIMULATION VALUES. 
2). VALUES IN PARENTHESES ARE VALUES FROM THE SMOOTHED 
DYNAMIC TRANSITION MAP. 
3). THE TWO VARIALBLE LINEARLY REGRESSED EQUATION FOR 
MATRIX COEFFICIENT Al2 IS: 


Al2 = -1. O*EXP(-2.8415*(NSL/NSB) + 
7.9978*(NGL/NGB) - 4. 4662) 
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TABLE 6 


A13 MATRIX COEFFICIENT CORRELATION DATA 


57 000 7RPM 


35000 


32000 


30000 


25000 


23000 


22000 


20000 


17000 


15000 


NOTES: 


NS = 500 RPM 1000 P00 2000 2500 
4410. 37 3869.54 pogo. 035 297 Cat 2613.44 
4024. 89 Ses 33 5093. 29 2718. 36 Zo8o. UL 
S08. 91 2073.62 2701. 10 2509 MOU 207926 

(4500, 00) (3050.00) (2200. 00) (1990.00) (2000. 00) 
SZ Oe 22 2809.54 2465.01 Zeca LOIS 5 
ero Oo od) 3020.00) "C2 P00; 00) (1920.00) (1900. 00) 
2547. 69 ZL 28 Pot i] 1720 oo 1509.68 
Cy o0200) C2550. 00) (1800.00) (1760.00) (1800.00) 
2325.02 209. 91 oO. 6 157 U2 vey? 3 
Des OO 1948. 72 1709.76 1500509 1316.14 
(1800.00) (1550.50) (1450.00) (1430.00) (2000. 00) 
2026296 17 Gino 9 IL S16) OS) P66. 96 ZO tet 
176) 11: 1550.41 leu. 29 1193.48 1047. 13 
NGrZ. 66 1414. 90 Po. 39° 1089. 17 ZS) Sasol 


VALUES WITHOUT PARENTHESES ARE SIMULATION VALUES. 

VALUES IN PARENTHESES ARE VALUES FROM THE SMOOTHED 
DYNAMIC TRANSITION MAP. 
THE TWO VARIALBLE LINEARLY REGRESSED EQUATION FOR 
MATRIX COEFFICIENT A113 IS: 


Al3 = -1. O*EXP(-0. 45788*(NSL/NSB) + 
1.1890*(NGL/NGB) + 6.8305) 
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TABLE 7 


A21 MATRIX COEFFICIENT CORRELATION DATA 


NS = 500 RPM 1000 1500 2000 2500 

NG = 37000 RPM eelel O29 Omr3 OmS0 0.45 
35000 O772 0.74 0556 0.45 O° 33 
32000 O567 O52 0359 On 29 0.22 


CO 0}) (0.50) (0. 40) (0730) (0.239) 


30000 O25: 0.40 0229 0.22 0.6 
CoO) (0.40) (0. 30) (O7Z6)) CO. Tea 


25000 Om 0.18 O315 0.09 OR Oe, 
C0230) (0,20) COm IG) (0, 10) (0.2oF 


23000 0.18 OMe On? 0. 08 0. us 


22000 0.14 O..09 0. 08 0. 08 O21 
COT10) COPED COMO) (0. 10) (O.GH 


20000 O7209 0.707 C207 OF 08 0. Ts 
17000 O05 0. 06 C209 Os 0.28 
15000 O.505 Ones On13 Onze 0.20 


NOTES: 1). VALUES WITHOUT PARENTHESES ARE SIMULATION VALUES. 
2). VALUES IN PARENTHESES ARE VALUES FROM THE SMOOTHED 
DYNAMIC TRANSITION MAP. 
3). THE TWO VARIALBLE LINEARLY REGRESSED EQUATION FOR 
MATRIX COEFFICIENT A21 IS: 


A21 = (-1.5312E-01)*((NSL/NSB)**2.0) + (-9.5315E-01)*(NSL/NSB)*(NGL/NGB) 
(1.5745)*CCNGL/NGB)**2.0) + (5.1810E-01)*(NSL/NSB) + 
(-1.6232)*(NGL/NGB) + 4.6015E-01 


TABLE 8 


A22 MATRIX COEFFICIENT CORRELATION DATA 


NS = 500 RPM 1000 1500 2000 2500 

NG = 37000 RPM -1.46 -1.44 aaa S -1.41 =i 9 
35000 aE, oO -1. 34 ale oL “l, 22 =i, 28 
32000 eile 2 1 aL 9 = =1..16 -1.14 


e105) a2 0) G12 25) Gal. 25) C120) 


30000 =1'10 = 0S ale 7 = O0 -1.04 
(lO, ES) ele 10:) Cor, 10) Colle op, 00) 


25000 SURG. =O 835 a0. G2 =0)230) = O.a7 2 
C- 055) C010) C=0°65) C0760) (057 O) 


23000 Swe 7D S027 3 SUT 2 =0°70 20706 


22000 -0. 69 -0. 68 -0. 67 -0.65 -0. 63 
(=-0.75) (-0.80) €-0.70) (-0.60) (-0.50) 


20000 Ooo SOL oO Ooo 7 =0\55 ORO 
17000 -0.45 -0. 43 -0.41 =O o0 =Va36 
15000 -0. 34 a0. 33 20.31 =O =o. 23 


NOTES: 1). VALUES WITHOUT PARENTHESES ARE SIMULATION VALUES. 
2). VALUES IN PARENTHESES ARE VALUES FROM THE SMOOTHED 
DYNAMIC TRANSITION MAP. 
3). THE TWO VARIALBLE LINEARLY REGRESSED EQUATION FOR 
Math Ex CORFE ICTENT =A22 15: 


A22 = (5.6875E-02)*(NSL/NSB) + (-1. 3166)*(NGL/NGB) 
rc. 3. 9002 E= Om 
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TABLE 9 


A23 MATRIX COEFFICIENT CORRELATION DATA 


NS = 500 RPM 1000 1500 2000 2500 
NG = 37000 RPM 55 220 | Ze OZ 3.41 4.43 
35000 Ze Zo 3205 4.72 6.14 
32000 325.1 4.56 Dee Toa 10. 04 


(4. 00) (4.00) (4. 00) COO U) (10.005 


30000 4.8/7 6. 33 S223 LO aval 13.33 
(4. 00) (4.00) CIO) (15.00) (153003 


25000 P03 14. 35 18. 66 24.27 31Sr 
(5.00) (20. 00) (24550) (26. 00) ( 2850s 


23000 Looe 19: 90 25-89 33368 43.79 


22000 Te-.02 30.49 59.00 51.256 67. 02 
(25750) (30700) (32700) (34. 00) (35. 0GF 


20000 25.00 Soe 42, 29 55 e00n 7 lS 
17000 40. 85 Doe OoralO 89366 116.538 
15000 56. 66 Toe 96-08 124. 69 162.9he 


NOTES: 1). VALUES WITHOUT PARENTHESES ARE SIMULATION VALUES. 
2). VALUES IN PARENTHESES ARE VALUES FROM THE SMOOTHED 
DYNAMIC TRANSITION MAP. 
3). THE TWO VARIALBLE LINEARLY REGRESSED EQUATION FOR 
MATRIX COEFFICIENT A23 IS: 


A23 = -1. O*EXP(0.92911%*(NSL/NSB) - 
4. 2459%(NGL/NGB) + 6.2290) 
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APPENDIX G 


SIMULATION VS. DATA RUNS 


This appendix contains the results of the three computer 
simulations used to validate the Smoothed Dynamic Transition 


Map. Figures 20-25 document the results for both NG and NS. 
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